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PROBLEM 


Provide  analysis  and  synthesis  of  Navy  design  problems.  Specifically, 
develop  a  capability  at  NELC  for  using  mathematical  programming  (MP)  as  an 
aid  to  engineering  design. 


RESULTS 

1.  The  overall  software  capabilities  at  NELC  for  numerically  solving 
various  types  of  MP  problems  -  initiated  or  developed  under  this  problem  - 
are  discussed  in  the  report  proper.  Applications  of  integer  programming  are 
given  in  Appendix  1.  User’s  guides  and  FORTRAN  codes  f'^r  solving  some 
classes  of  MP  problems  are  given  in  Appendix  2. 

2.  Practical  guidelines  are  given  for  applying  MP  methods. 


RECOMMENDATIONS 

1 .  Maintain  a  continued  effort  to  keep  the  mathematical  programming 
software  current.  Monitor  research  literature  for  new  developments  in  non¬ 
linear  programming  and  integer  programming. 

2.  Conduct  ongoing  seminars  or  in-house  classes  to  inform  practicing 
scientists  and  engineers  of  the  utility  of  MP. 

3.  Review  Navy  engineering  design  problems  for  possible  application 
of  MP  techniques. 
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INTRODUCTION;  SCOPE  OF  REPORT 


There  are  several  problems  involved  in  the  development  of  faithful 
mathema  dcal  models  of  real-world  processes.  The  processes  are  in  general 
nonline  rr.  The  modeling  equations  are  frequently  incomplete.  Conditions 
are  known  only  within  limits.  Often  the  best  approach  to  these  problems  is 
via  mathematical  programming  (MP).^‘^ 

MP  is  a  distinct  discipline  -  it  exists  independently  of  computer  pro¬ 
gramming.  Prior  to  the  age  of  the  high-speed  computer,  however,  some  of 
the  original  algorithms  for  the  solution  of  MP  problems  were  too  cumbersome 
to  be  of  real  use.  The  advent  of  the  modern  digital  computer  hgs  made  the 
solution  of  many  types  of  MP  problems  feasible  and  has  stimulated  the  search 
for  better  algorithms. 

This  report  is  chiefly  concerned  with  solving  MP  problems  -  with  the 
computational  stage  of  MP.  It  is  intended  as  a  guide  for  the  usage  and  appli¬ 
cation  of  available  MP  c^  mputer  codes. 

THE  MATHEMATICAL  PROGRAMMING  PROBLEM  defines  the 
general  MP  problem. 

MATHEMATICAL  PROGRAMMING  CAPABILITIES  AT  NELC  lists 
and  evaluates  the  MP  codes  operational  on  the  NELC  IBM  360/65  computer, 
and  will  be  of  interest  to  the  user  who  has  an  MP  problem  in  final  form,  ready 
to  solve.  He  can  choose  the  appropriate  code  from  the  list  and  obtain  the 
card  deck  and  user’s  guide  from  the  NELC  program  library.  Computer  Sciences 
Department. 

DESIGNING  WITH  MAIHEMATICAL  PROGRAMMING  AS  AN  AID 
provides  guidelines  for  modifying  the  MP  problem  when,  m  its  first  form,  it 
is  cumbersome;  or  when  there  is  not  information  enough  to  start  computation 
(for  example,  an  initial  feasible  point  is  lacking);  or  when  the  available  codes 
do  not  yield  all  the  needed  information  (for  example,  postoptimal  analysis). 


;.  See  APPENDIX  3:  REFERENCES. 
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THE  MATHEMATICAL  PROGRAMMING  PROBLEM 


The  basic  problem  of  MP  is  to  develop  an  algorithm  for  finding  the 
minimum  of  a  scalar-valued  function  of  n  real  variables  that  satisfies  a  set  of 
auxiliary  conditions  called  constraints.  Stated  in  mathematical  terms,  the 
problem  becomes; 

Let  f(X] - ,  x„)andgj(xj, . . . ,  x„), . . .  ,gn,(xi, . . . ,  x„)  be 

scalar-valued  functions  of  the  n  variables  Xj, . . . ,  Xj^.  Then  we  wish  to  find 
variables  which 

minimize  f(xj,  X2, . . . ,  Xj^) 
subject  to 

gj(xi,X2,  ...,x„)>0  (1) 


gm(xi,X2,  ...,Xn)>0 

The  above  problem  is  known  variously  as  the  ‘general  mathematical  program¬ 
ming  problem,’  the  ‘constrained  optimization  problem,’  and  the  ‘nonlinear  pro¬ 
gramming  problem.’  For  the  sake  of  convenience,  we  call  it  problem  (1). 

In  problem  (1 ),  f  is  called  the  obje«;tive  or  cost  function  and  the  gj  are 
called  the  constraints.  We  also  refer  to  f  and  gj  as  the  problem  functions. 

We  denote  the  vector  (x  j ,  X2,  .  • x„)  by  x^  (where  T  denotes  the  matrix 
transpose)  and  call  the  set  of  all  vectors  x  which  satisfy  gj(x)  >0  for  all  i= 1 , 

. . . ,  m,  the  constraint  set  or  feasible  set.  The  problem  is  said  to  be  consistently 
posed  if  the  constraint  set  is  nonempty.  We  note  that  finding  the  maximum  of 
a  function  f  is  equivalent  to  finding  the  minimum  of  -f. 

Consider  the  following  example  (fig.  1 ): 

Minimize  f(xj,  X2)  =  (Xj-2)^  +  <X2-1)^ 
subject  to 

gj(Xj,X2)  =  X2-x^>0  (2) 


g2(X|,X2)  =  -xj  -X2  +  2^0 


Since  f(xj ,  X2)  is  the  sum  of  squares,  the  minimum  occurs  at  xj  =  2,  X2  =  1 
However,  the  point  (2,1 )  is  not  in  the  constraint  set  defined  by  g|  and  g2. 
The  obvious  (in  this  example)  candidate  is  the  point  (1,1),  which  in  this 
straightforward  problem  is  the  constrained  minimum.  As  m  and  n  get  larger, 
the  problem  becomes  significantly  more  difficult  to  solve. 
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Figure  1.  Constraint  set:  MP  problem. 


We  now  turn  to  a  discussion  of  computer  codes  which  can  solve  MP 
problems,  for  different  classes  of  problem  functions. 

MATHEMATICAl-  PROGRAMMING  CAP  ABILITIES  AT  NELC 

UNCONSTRAINED  PROBLEMS 

The  unconstrained  MP  problem  is  stated  as: 

Minimize  f(xj,  X2, . . . , 

where  f  is  a  scalar-valued  function  of  the  n  variables  (xj,  X7, .  .  .  ,  Xj^)^  =  x. 

Gradient  methods  and  direct  search  techniques  are  the  two  basic 
approaches  to  numerically  solving  the  unconstrained  problem.  The  gradient 
of  f  at  X,  denoted  by  Vf(x),  is  defined  to  be  the  following  vector 

Vf(x)  =  (3f(x)/ax|, . . .  ,  3f(x)/9x„)^ 

Gradient  methods  use,  in  some  way,  the  following  facts; 

1.  At  the  minimum  x*  of  f,  we  have  Vf(x*)  =  0. 

2.  If  Vf(x)  ^  0,  then  -Vf(x)  points  in  the  direction  of  steepest  descent. 


I 


tA  finer  classification  would  be  direct  methods,  gradient  methods,  and  those  methods 
involving  the  matrix  of  second  derivatives.  We  feel  that  the  first  two  are  the  most  useful 
for  applications. 
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This  analytic  information  makes  the  gradient  methods  fast  and  enables  the 
computer  codes  to  compute  meaningful  error  information.  To  use  gradient 
methods,  we  must  have  the  gradient  of  f  available  analytically  or  have  a  numeri 
cal  way  to  compute  it.  Direct  search  methods  eliminate  this  need  for  the 
gradient  and  rely  only  on  the  behavior  of  the  objective  function  in  seeking  out 
the  minimum.  Typically,  direct  methods  evaluate  the  cost  function  many 
more  times  than  gradient  methods  in  minimizing  the  sane  test  function.  In 
minimizing  the  Rosenbrock  test  function  (see  Appendix  2),  the  direct  search 
routine  ZANGWL  requires  325  function  evaluations,  while  the  gradient 
method  CONJGT  requires  only  71  combined  function  and  gradient  evalua¬ 
tions,  in  finding  the  optimum  to  within  the  same  accuracy.  This  is  a  trade¬ 
off  a  user  must  make  if  he  can  choose  between  a  gradient  method  and  a  direct 
method. 

Direct  methods  do  not  rest  on  so  firm  a  mathematical  foundation  as 
gradient  methods  do,  and  most  direct  methods  are  proved  to  converge  for 
only  special  functions.  However,  they  have  been  useful  m  practice,  since  the 
objective  function  in  many  applications  is  complicated  or  its  gradient  is  not 
available.  It  is  generally  simpler  to  code  a  problem  for  a  direct  method,  which 
allows  for  faster  implementation  on  the  computer. 

We  present  the  following  example  of  posing  a  two-point  boundary 
value  problem  (TPBVP)  as  an  unconstrained  MP  problem,  to  illustrate  both 
an  application  of  MP  and  the  need  for  good  direct  search  methods.  The  prob¬ 
lem  is  to  find  an  n-dimensional  vector  function  y(t)  which  satisfies 

y  =  h(t,y)  (a<t<b) 

with 

yi(a)^  qia  (i=l,2, . . .  ,j<n) 

and 

yi(b)  =  qib  (i=j+l,...,n) 

In  general,  this  problem  has  no  closed-form  solution,  and  in  some  cases  no 
solution  at  all.  However,  since  it  frequently  arises  in  applications,  either  a 
numeiica!  estimate  of  the  solution  is  desired  or,  if  no  solution  exists,  a  function 
which  comes  close  (in  some  sense)  to  solving  the  problem  is  desired.  With  this 
in  mind  we  pose  the  following  MP  problem.^ 

/  n 

Minimize  fix  j,X2,...,Xn_j)=  I  ^  (yi(b)-qjj,)^  |  (3) 

_ / 

^Rosen^  discusses  the  same  problem  and  obtains  approximate  solutions  using  linear  pro¬ 
gramming  techniques.  His  approach  requires  a  great  deal  of  equation  manipulation  before 
the  linear  programming  techniques  can  be  applied.  Unfortur;ately,  no  comparison  between 
the  two  methods  has  been  made. 
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where  the  numbers  yj(b),  i=j+ 1, . . . ,  n  are  numerically  computed  as  follows. 
For  a  given  (Xj, . . .  ,  solve  the  following  initial-value  problem  over  the 
interval  [a,b] . 

y  =  h(t,y)  (a  <  t  <  b) 

where 

y(a)  =  (qia,  ...,qja,xi,...,x„.j) 

the  last  n-j  components  of  the  solution  obtained  at  t=b  of  this  problem  are 
used  in  the  objective  function  for  yj(b),  i=j+ 1, .  .  .  ,  n. 

In  this  problem  there  is  no  analytic  expression  for  the  objective 
function  f,  from  which  Vf  can  be  derived.  Thus,  to  numerically  solve  this 
unconstrained  problem  (3),  either  a  direct  search  method  must  be  used  or 
Vf  must  be  numerically  calculated  via  a  differencing  routine.  We  recommend 
the  former  as  reliable  and  easy  to  use,  and  discuss  some  of  the  drawbacks  of 
the  latter  in  DESIGNING  WITH  MATHEMATICAL  PROGRAMMING  AS 
AN  AID,  r radient  approximation.  Dejka^  discusses  a  similar  TPBVP  and 
uses  a  direct  search  method  to  solve  a  related  MP  problem. 

We  briefly  discuss  computer  routines  available  from  the  NELC  pro¬ 
gram  library  for  solving  the  unconstrained  p'-oblem.  The  user’s  guides  for 
these  routines  provide  ample  background  information  ana  details  for 
implementation. 

Gradient  methods  from  the  library  are  FP,  CONJGT,  SOREN,  FMFP, 
and  FMCG.  FP  and  CNJGAT  originally  were  programmed  and  used  by 

O 

Winterbauer  to  solve  a  parameter  selection  problem  for  a  sonar  signal  equa¬ 
tion,  but  they  are  general-purpose,  unconstrained,  MP  codes.  FP  and  CNJGAT 
are  based  on  the  methods  of  Fletcher  and  Powell^  and  Fletcher  and  Reeves,'^ 
respectively.  SOREN  is  a  modification  of  CNJGAT  which  has  converged 
faster  for  some  test  functions.  FMFP  and  FMCG  are  available  from  the  IBM 
Scientific  Subroutine  Package;*  *  we  have  not  tested  these  last  two  routines 
extensively  or  compared  them  with  the  other  gradient  methods. 

ZANGWL,  direct,  and  UNIV.AR  are  direct  search  codes  in  the 
library.  ZANGWL  is  discussed  in  Appendix  2,  and  DIRECT  and  UNI  VAR  are 
presented  in  reference  12.  These  three  routines  are  based  on  methods  presented 
in  references  13-15,  respectively,  and  were  programmed  at  NELC.  ZANGWL 
has  a  mat^'ematical  basis  for  convergence  similar  to  that  of  CNJGAT,  and  of 
the  three  direct  methods,  it  is  the  most  efficient  in  terms  of  the  1  >tal  number 
of  function  evaluations  lequired  to  minimize  a  function.  In  minimi/.ing  the 
Rosenbrock  function,  to  the  same  accuracy  and  from  the  same  initial  point, 
the  number  of  function  evaluations  were  ZANGWLt525),  D1RFCT(705),  and 
UNIVAR(2303).  The  extraordinary  number  of  function  evaluations  clearly 
makes  UNIVAR  unacceptable,  but  ZANGWL  should  not  be  selected  over 
DIRECT.  DIRECT  makes  intermediate  searching  moves  in  a  much  more 


cautious  manner  than  ZANGWL,  which  makes  it  better  for  some  applications. 
This  is  discussed  in  DESIGNING  WITH  MATHEMATICAL  PROGRAMMING 
AS  AN  AID,  SUMT  and  Constraint  Transformation. 

Before  moving  onto  the  constrained  problem,  we  give  a  word  of 
caution  on  all  minimization  routines.  Each  method,  be  it  a  constrained  or  an 
unconstrained  code,  is  capable  of  finding  only  a  local  minimum  and  not  a 
global  minimum.  We  define  local  and  global  minima  in  the  next  lines.  Let 
the  objective  function  f(x)  be  defined  on  a  set  G  in  an  n-dimcnsional  vector 
space,  denoted  by  E".  Then  we  say  that  f  has  a  global  minimum  at  x*  (in  G) 
il  f(x*)  <  f(x)  for  all  x  in  G.  Note  that  we  do  not  exclude  G  =  E”.  x  in  G  is 
called  a  local  minimum,  if  for  all  x  sufficiently  close  (with  respect  to  some 
norm)  to  x,  and  also  in  G,  we  have  f(x)  <  f(x). 

Local  minima  can  occur  in  the  gradient  methods  because  the  condition 
that  Vf(x*)  =  0  is  only  necessary  and  not  sufficient  for  a  global  minimum.  In 
direct  methods,  only  local  information  about  the  surface  defined  by  the  objec¬ 
tive  function  is  available  to  the  routine.  This  characteristic  makes  direct 
methods  susceptible  to  stopping  at  a  local  minimum.  For  reasonable  certainty 
that  a  global  optimum  has  been  reached,  it  is  wise  to  restart  the  problem  from 
different  initial  points.  In  many  applications,  if  a  local  minimum  gives  a 
satisfactory  value  of  the  objective  function,  no  further  processing  is  necessary. 

Work  is  continuing  in  the  area  of  unconstrained  minimization  algorithms, 
with  refinements  to  the  above  methods  and  new  methods  appearing  regularly 
in  the  literature.  The  most  fruitful  and  accessible  sources  of  articles  on  the 
subject  aie  The  Computer  Journal,  Communications  of  the  Association  of 
Computing  Machinery,  SIAM  Review,  and  the  SIAM  journals  on  control, 
numerical  analysis,  and  applied  mathematics.  The  above  sources,  together  with 
Management  Science  and  Operations  Research,  contain  many  articles  on  the 
constrained  problem. 

CONSTRAINED  PROBLEMS 

We  return  to  the  discussion  of  problem  (1)  for  various  classes  of 
problem  functions.  The  following  types  of  mathematical  programming  prob¬ 
lems  are  discussed:  linear,  quadratic,  convex,  nonlinear  and  nonconvex,  and 
integer.  The  methods  for  solving  these  problems  make  explicit  use  of  the 
propeiii  o  of  the  problem  functions. 

When  an  unconst.ained  problem  is  solved,  the  codes  require  only  an 
initial  point  for  which  the  objective  function  is  defined.  Th's  requirement  is 
more  demanding  in  the  constrained  problem.  Depending  on  the  type  of  prob¬ 
lem  under  consideration,  the  user  can  be  required  to  provide  an  initial  feasible 
point,  as  a  starting  point  for  the  computation.  In  many  applications  an  initial 
feasible  point  is  Known  from  the  engineering  knowledge  of  the  problem.  How¬ 
ever,  if  the  constraints  are  numerous  or  crmplicated,  such  a  point  will  not  be 


obvious  and  a  preliminary  step  must  be  taken  prior  to  solving  the  problem. 

A  method  for  obtaining  an  initial  feasible  point  is  treated  in  DESIGNING  WITF 
MATHEMATICAL  PROGRAMMING  AS  AN  AID,  Initial  Points  and  Scaling. 
so  we  assume  that  one  is  at  hand  in  the  following  discussion. 

Each  class  of  constrained  MP  problem  is  described,  together  with 
computer  programs  which  can  solve  it.  Since  each  routine  to  be  discussed  has 
an  associated  user’s  guide,  we  confine  our  remarks  to  the  following  points: 

1 .  Is  an  initial  feasible  point  required? 

2.  Does  the  code  find  a  global  optimum? 

3.  Can  information  be  saved  for  possible  restarts  or  postoptimal 
analysis? 

4.  Is  the  routine  easy  to  use? 

5.  What  error  messages  are  given  if  tne  routine  fails  to  converge? 

LINEAR  PROGRAMMING  (LP) 

In  the  standard  linear  programming  problem  all  the  functions  in 
question  are  linear  and  the  problem  variables  are  constrained  to  be  non¬ 
negative.  We  write: 

Minimize  f(x)  =  CjXj  +  C2X2  +  + 

subject  to 

gl(x)  =  aiiXi+a,2X2  +  ...  +  a,nX„<b, 


"  ^ml'^l  +  ^in2’‘2  +  •  •  •  +  ^mn’^n  ^ 

This  LP  problem  can  be  written  in  matrix  notation  as: 

Minimize  z  =  c^x 
subject  to 

Ax<b 

where  the  <=  means  that  the  conesponding  components  of  the  vectors  are 
“less  than  or  equal  to.”  If  the  constraints  are  consistent,  then  the  simplex 
method  of  linear  programming  guarantees  that  a  global  optimum^  can  be 
found  in  a  finite  number  of  steps.  The  simplex  method  is  an  iterative  proce¬ 
dure  and  generates  “basic  feasible  solutions”  at  each  iteration,  which  decrease 
z.  To  produce  these  solutions,  a  “basis  inverse”  matrix  is  calculated.  The 

^An  unbounded  solution  can  also  be  detected  in  a  finite  number  of  steps. 


preceding  brief  comments  serve  only  to  associate  the  terms  “basic  feasible 
solution”  and  “basis  inverse”  witi\  the  simplex  method;  references  16  and  17 
treat  the  simplex  method. 

Tlie  most  complete  code  for  numerically  solving  the  LP  problem  is  the 
IBM  Mathematical  Programming  System^®  (MPS/360).  MPS/360  is  based  on 
a  modification  of  the  simplex  method  and  will  either  solve  the  LP  problem  or 
indicate  that  no  solution  exists.  This  code  does  not  require  that  the  problem 
variables  be  nonnegative,  and  treats  upper  and  lower  bounds  on  the  variables 
as  special  constraints.  Separable  programming  problems  (a  special  nonlinear 
MP  problem)  can  be  solved  with  this  routine.  MPS/360  is  capable  of  solving 
problems  of  up  to  4095  constraints  and  “virtually  an  unlimited  number  of 
columns.”*  ^  It  is  currently  stored  on  disk  pack  NELC05  at  the  NELC 
Computer  Center. 

No  initial  point  is  required  to  begin  the  computation;  however,  the 
option  exists  to  start  the  problem  from  a  user-supplied  basis  inverse.  MPS/360 
has  its  own  control  language,  which  provides  a  variety  of  capabilities.  The 
user  is  afforded  several  postoptimal  analysis  procedures  and  can  access  the 
current  basis  inverse  for  future  restarts.  This  control  language  is  straight¬ 
forward  to  use  and  provides  some  looping  and  branching  capability.  A  variety 
of  messages  are  output  to  the  user  in  the  course  of  computation.  They  are 
fully  explained  in  the  message  manual.*^ 

The  chief  drawback  of  this  program  is  the  format  of  the  input  data. 

It  requires  each  element  of  the  arrays  c,b,  and  A  to  have  a  “row  name”  and  a 
“column  name”  for  identification.  This  has  proved  cumbersome  for  scientific 
and  engineering  work.  A  FORTRAN  program,  DATAPREP,  is  available  to 
reduce  the  data  arranged  in  compact  matrix  notation  to  a  format  acceptable 
to  MPS/360. 

In  many  applications  a  lineai  programming  problem  must  be  solved 
repeatedly  as  part  of  a  larger  problem.  The  READCOMM^®  facility  of  MPS/360 
allows  the  main  program  to  be  used  in  an  iterative  fashion  as  a  subroutine. 
READC'OMM  enables  the  user  to  supplement  the  standard  control  language 
with  FORTRAN  procedures;  for  example,  DATA  PREP.  Rosen^  and  Griffith 
and  Stewart^  *  have  examples  of  using  a  linear  programming  code  in  an  itera¬ 
tive  way. 

Previous  large-scale,  efficient  LP  codes  were  geared  to  commercial 
applications  and  required  a  great  deal  of  modification  for  efficient  scientific 
and  engineering  use.  The  READCOMM  facility  has  made  a  powerful  program 
easily  available  for  a  wide  range  of  specialized  applications. 

QUADRATIC  PROGRAMMING  (QP) 

This  type  of  problem  is  the  next  order  of  difficulty.  A  quadratic 
cost  function  is  minimized  subject  to  linear  constraints: 


Minimize  f(xj,  X2, . . . ,  x^) 
subject  to 


a,iXi+ai2X2  +  .-  +  ainXn<b, 


(4) 


aml^l+W2+-  - 

where  f  has  one  of  two  forms  - 

f(xj,X2, . . . ,  Xjj)  =  cTx  + x^Bx  (5) 

or 

f(xj,  X2, . . . ,  Xj^)  =  llHx -ell  (6) 

B  and  H  are  n-by-n  and  k-by-n  matrices,  respectively,  and  c  and  e  are  n-,  and 
k-dimensional  vectors,  respectively.  The  norm  of  a  vector  y,  denoted  by  lly  II, 
is  given  by 


At  present  only  the  minimum  norm  problem  (equation  6)  can  be 
solved  at  Nij-LC.  The  program  which  does  this  is  QPH  AN  SON.  This  routine 
was  written  by  R.  J.  Hanson,  of  the  Jet  Propulsion  Laboratory,  Pasadena,  and 
uses  a  numerically  stable^  version  of  Rosen’s^^  gradient  projection  algorithm. 
The  method  guarantees  that  a  global  Optimum  will  be  found  for  a  consistent 
problem.  The  routine  is  reported  to  have  worked  well  on  examples  from  the 
areas  of  curve  fitting  and  approximation  of  solutions  to  linear  integral 
equations. 

The  routine  is  described  in  reference  23,  and  an  NELC  user’s  guide  is  in 
preparation.  Tlie  code  is  operational,  but  it  is  still  unpolished  in  respect  to 
user-oriented  input  and  output.  QPHANSON  does  not  require  an  initial 
feasible  poini  nor  does  it  have  an  option  to  accept  a  good  approximation  of 
the  optimum.  The  code  treats  equality  constraints  and  inequality  constraints 
separately.  It  can  solve  problems  of  up  to  60  constraints  (equality  and  in¬ 
equality  combined)  and  30  variables.  No  experiments  have  been  done  to 
determine  whether  this  is  a  hard  and  fast  upper  bound  on  the  problem  size. 

The  program  suffers,  from  lack  of  good  error  and  timing  messages.  If 
the  routine  fails  to  converge,  no  messages  are  given  as  to  the  possible  cause. 
Also,  no  provisions  are  made  to  identify  inconsistent  quadratic  programs. 

The  user  is  on  his  own  with  QPHANSON. 


^An  algorithm  0  is  numerically  stable  if  the  errors  in  the  input  data  are  approximately 
equal  to  the  round-off  errors  generated  by  the  computations  of  0, 


If  a  QP  problem  occurs  with  equation  (5)  as  the  cost  function,  and  the 
matrix  B  is  positive-definite  or  positive-semidefinite,^^  then  Hanson^^  presents 
a  method  for  transforming  this  QP  problem  into  a  minimum-norm  QP  prob¬ 
lem.  This  minimum-norm  problem  is  solved  with  QPHANSON  and  then  an 
inverse  transformation  is  made.  Presently  this  must  be  done  by  the  user. 
QPHANSON  is  currently  being  modified  to  make  this  transformation 
automatically. 

Codes  to  solve  quadratic  programs  are  not  as  well  polished  or  as  highly 
developed  as  those  for  LP.  Unless  the  demand  increases  for  good  QP  routines, 
the  user  will  have  to  write  his  own  code  or  be  content  wit  i  the  experimental 
models. 


CONVEX  PROGRAMMING 

Before  turning  to  the  convex  programmirg  problem,  we  make  some 
preliminary  definitions. 

CONVEX  SET.  A  set  G  in  E”  is  said  to  be  convex  (fig.  2)  if  for  any 
two  points  Xj  and  X2  contained  in  G  we  have  Xxj  +  (1  -\)x2  contained  in  G, 
for  all  X  in  (0,1). 


nonconvex  sets 


Figure  2.  Examples  of  convex  and  nonconvex  sets. 


symmetric  matrix  B  is  positive -definite  (semidefinite)  if  for  every  x  #  0  we  have 
xT'Bx  >0,(^0). 
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CONVEX  FUNCTION.  A  scalar-valued  function  f  defined  on  a  convex 
set  G  in  E"  is  said  to  be  convex  if  for  any  two  points  X|  and  X2  in  G 

f(XX|  +  (1  -X)x2)  <  Xf(xj)  +  (1  -X)f(x2) 
for  all  X  in  (0,1 ). 

Linear  functions,  and  the  quadratic  cost  function  (equation  6)  with  B 
positive-semidefinite,  are  examples  of  convex  functions.  A  theorem  of  interest 
states  that  if  the  constraint  set  of  an  MP  problem  is  defined  by  convex  functions, 
then  it,  too,  is  convex.  More  precisely,  if  g|(x), . . . ,  are  convex  func¬ 
tions,  then  the  set  of  all  x  for  which  gj(x)  <  0, . .  . ,  g^^fx)  <  0  holds  simul¬ 
taneously  is  a  convex  set.  The  constraint  sets  of  linear  and  quadratic  programs 
are  convex. 

With  these  facts  in  hand  we  state  the  convex  programming  problem. 

Minimize  f(x) 
subject  to 

gj(x)<0  (i=l,...,m) 

where  the  functions  t  and  gj  are  convex  functions  of  x  =  (xj , . .  .  ,  x^^)* . 

. .  great  deal  of  work  has  been  done  with  convex  programming  and 
the  theory^^  can  guarantee  convergence  for  some  computational  methods. 

Each  method  is  valid  for  specific  requirements  on  the  problem  functions.  In 
this  section  we  assume  that  the  gradients  of  all  functions  exist  and  are  con¬ 
tinuous.  The  central  problem  in  solving  convex  programs  is  not  so  much 
theoretical  difficulty  but  rather  the  obtaining  of  rapid  convergence  of  numeri¬ 
cal  schemes.  Even  though  some  QP  problems  are  convex,  it  may  be  more 
efficient  to  use  a  routine  tike  QPHANSON  to  solve  them  rather  than  treat 
them  as  convex  programs.  Another  practical  difficulty  with  convex  program¬ 
ming  is  the  identification  of  convex  functions.  If  the  function  has  a  complicated 
analytic  expression,  it  can  be  difficult  to  classify  it  as  convex.  The  methods 
for  solving  convex  programs  will  not  in  general  completely  hang-up  if  the  data 
are  not  convex,  but  the  significance  of  s  ;ch  results  should  be  judged  in  terms 
of  the  user’s  problem  formulation. 

The  first  library  routine  which  solves  the  convex  programming  prob¬ 
lem  is  the  subroutine  CONVEX,  which  was  developed  by  Hartley  and 
Hocking^^  at  Texas  A&M.  The  routine  makes  a  linear  approximation  to  the 
functions  in  question  and  then  uses  a  simplex-like  procedure  to  move  to  the 
optimum.  In  making  the  linear  approximations,  the  routine  requires  a  user- 
supplied  subroutine  which  computes  the  gradients  of  the  cost  function  and 
the  nonlinear  constraints. 

CONVEX  does  not  require  an  initial  feasible  point:  however,  the 
option  does  exist  to  start  from  a  given  point.  In  addition.  CONVEX  produces 


a  current  feasible  point  and  a  basis  inverse  at  the  end  of  each  i  ^ration  for 
possible  restarts.  The  format  of  the  input  data  is  straightforward  and  suitable 
to  scientific  and  engineering  work;  however,  care  should  be  exercised  in  the 
organization  of  the  data  for  any  upper  and  lower  bounds  on  the  problem 
variables.  CONVEX  requires  that  the  constraint  data  be  input  in  three 
groups  -  the  upper  and  lower  bounds  on  the  variables,  the  linear  inequalities, 
and  the  nonlinear  convex  constraints.  This  feature  makes  it  possible  to  con¬ 
veniently  solve  quadratic  programs  with  convex  cost  functions.  No  compari¬ 
sons  between  QPHANSON  and  CONVEX  have  been  made  on  solving  quadratic 
programming  problems. 

CONVEX  suffers  from  the  lack  of  good  error  messages  and  analysis 
in  the  event  of  an  inconsistent  problem  or  any  numerical  difficulties.  No 
investigation  of  the  numerical  stability  of  the  method  or  of  timing  or  accuracy 
benchmarks  for  large  problems  has  been  reported.  A  convex  problem  with 
60  constraints  and  60  variables  is  the  largest  which  can  be  solved  without 
program  modifications.  Because  the  linear  constraints  are  treated  separately, 
it  is  likely  that  larger  problems  can  be  solved  if  the  number  of  nonlinear  con¬ 
straints  is  not  too  great.  Future  work  should  investigate  this  possibility. 

The  second  routine  for  solving  the  convex  programming  problem  is 
Experimental  SUMT.  This  method  is  theoretically  convergent  for  convex 
data,  but,  since  it  also  has  provisions  for  nonconvex  programs,  we  postpone 
discussion  of  it  until  the  next  section. 

There  is  a  special  subclass  of  convex  programs  for  which  a  global 
optimum  can  be  found  with  the  linear  programming  code  MPS/360.  These 
are  separable  programming  problems,  which  are  defined  as  follows: 
n 

Minimize  z  =  ^  fj(Xj)  (7) 

j=l 


subject  to 


JgijfXjXbi 

j=I 


(i=l,...,m) 


Note  that  the  objective  function  and  the  constraints  are  sums  of 
functions  of  the  single  variables  Xj;  that  is,  there  are  no  “cross  product”  terms. 
This  allows  each  nonlinear  function  to  be  replaced  by  a  polygonal  approxima¬ 
tion,  and  reduces  problem  (7)  to  a  form  which  can  be  solved  by  MPS/360.  The 
MPS/360  user’s  manual  gives  the  appropriate  details  and  examples  of 
solving  separable  programming  problems.  If  the  separable  functions  are 
convex,  then  MPS/360  will  find  a  global  optimum  to  the  approximation 
problem.  The  use  of  successive  approximations  causes  the  global  optima  of 
the  approximation  problems  to  converge  to  the  optimum  of  the  original 
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problem.  The  method  can  tolerate  some  nonconvex  functions  but  may  stop 
at  a  local  optimum. 

This  technic  ie  can  also  be  used  for  separable  convex  programs  too 
large  for  CONVEX.  (CONVEX  may  be  more  efficient  if  the  problem  is  not 
too  large  -  unfortunately,  no  experimental  evidence  is  available  to  aid  in 
making  the  selection.) 

NONLINEAR,  NONCONVEX  PROGRAMMING 


This  last  class  of  MP  problems  is  composed  of  all  the  problems  which 
are  not  necessarily  linear  or  necessarily  convex.  In  the  statement  of  problem 
(1 ),  no  requirements  were  made  on  the  functions  in  question  other  than  the 
assumptions  that  the  objective  function  would  be  defined  for  all  feasible  x, 
and  the  gj  would  be  defined  for  This  statem  ;nt  of  problem  (1 )  is  much 
too  general  to  be  of  use.  To  have  any  hope  of  obtaining  a  solution,  we  must 
put  some  restrictions  on  the  problem  functions.  The  three  computer  codes 
which  we  discuss  require  that  the  gradients  of  all  the  functions  in  problem  ( 1 ) 
exist  and  be  continuous.  Although  these  conditions  are  stringent  from  a 
mathematical  point  of  view,^  they  do  not  provide  a  base  for  ar.  MP  algorithm 
The  following  example  illustrates  one  of  the  difficulties  of  nonlinear, 
nonconvex  programming.  Since  the  problem  functio.ns  are  possibly  nonconvex, 
complicated  constraint  sets  can  be  generated  (fig.  3). 


(x)  =  0 


9^(x)  =  0 


Figure  3.  Disconnected  constraint  set. 

In  this  example,  if  the  initial  feasible  point  x®  is  in  one  component  of  the 
constraint  set  and  the  optimum  x*  in  another  component,  then  the  routine 
cannot  move  from  x®  to  x*  while  keeping  intermediate  points  feasible.  Thus, 
the  most  that  general  MP  computer  routines  guarantee  is  a  local  minimum. 

One  of  the  most  sophisticated  routines  at  NELC  is  the  Ricochet 
Gradient'^  method.  This  is  an  IBM  SHARE  routine  which  requires  that  the 

^If  the  constraint  set  is  bounded,  then  these  conditions  are  sufficient  for  a  global  optimum 
to  problem  (1)  to  exist.  This  existence  theorem  gives  no  method  for  finding  the  optimum, 
which  can  be  difficult  in  the  general  problem. 


gradients  of  both  the  objective  function  and  the  constraints  exist  and  be  con¬ 
tinuous.  Tlie  method  requires  an  initial  feasible  point  and  begins  by  moving 
down  tlie  gradient  of  the  objective  function  until  a  constraint  is  reached.  The 
program  “ricochets”  and  traverses  on  the  objective  function  surface  across 
the  feasibility  region  to  the  opposite  constraint.  A  triangle  is  then  constructed 
with  this  traverse  line  as  its  base  and  its  apex  in  the  direction  of  the  gradient 
of  f.  The  next  step  is  made  along  the  line  from  the  base  to  the  apex.  The 
method  terminates  either  on  a  small  step  size  or  when  no  ricochet  is  possible. 

The  user  must  supply  codes  to  calculate  the  cost  function  and  the 
constraints  and  their  gradients.  An  initial  feasible  point  must  also  be  provided. 
The  program  has  no  options  for  restarts  or  postoptimal  analysis;  in  such  cases 
the  problem  is  simply  rerun  -  the  known  best  point  is  used  as  initial  data. 

This  routine  is  capable  of  producing  a  tremendous  amount  of  output 
information.  Once  the  method  for  controlling  this  output  has  been  mastered, 
the  user  has  access  to  a  variety  of  information,  which  can  be  of  great  use  in 
solving  a  nonlinear  problem.  The  accompanying  user’s  guide"^”  provides 
detailed  documentation  on  the  code  and  the  underlying  method.  A  supple¬ 
mental  user’s  guide  (Appendix  2)  reports  the  results  of  some  test  examples 
and  gives  a  sample  deck  setup  for  output  control.  This  program  has  proved 
reliable  and,  after  a  bit  of  experience,  easy  to  use. 

A  second  computer  code  for  the  general  problem  is  NELC  FESDIR. 

This  program  is  not  as  sophisticated  as  the  Ricochet  Gradient  routine,  but  it 
can  easily  be  modified  for  special  applications.  The  user’s  guide,  complete 
with  results  on  test  examples,  appears  in  Appendix  2. 

The  final  code  available  for  solving  the  nonhnear  problem  is  Experi¬ 
mental  SUMT,  which  was  written  at  Research  Analysis  Corporation  by 
G.  P.  McCormick,  et  a!.^^  SUMT  is  not  a  production  code  and  is  primarily 
used  as  a  research  tool  in  MP.  The  experimental  nature  does  not  lessen  its 
accuracy,  but  only  its  efficiency  and  speed.  The  code  has  modular  structure, 
which  allows  for  easy  user  modification  and  adaptation.  The  user’s  guide  is 
complete  with  test  examples,  although  some  handwritten  conections  and 
deletions  are  not  too  clear.  The  user  has  the  option  of  providing  an  initial 
point  himself  or  allowing  the  program  to  find  one.  The  option  also  exists  of 
having  SUMT  compute  the  gradients  by  a  differencing  scheme;  the  code  can 
also  check  user-supplied  gradients  for  errors.  SUMT  provides  timing  informa¬ 
tion  and  allows  for  user-controlled  output.  However,  the  output  can  be  con¬ 
fusing,  with  the  values  of  different  variables  appearing  under  the  same  headings. 
The  only  error  message  other  than  incorrectly  entered  data  is  a  warning  during 
computation  that  certain  estimates  indicated  the  problem  functions  to  be 
not  convex.  The  theoretical  background  for  SUMT  is  described  in  Fiacce 
and  McCormick. 

In  DESIGNING  WITH  MATHEMATICAL  PRCKiRAMMING  AS  AN 
AID,  SUMT  and  Constraint  Transformation,  we  discuss  methods  for  transforming 


a  constrained  problem  into  a  sequence  of  unconstrained  problems  (also  called 
SUMT).  This  requires  more  work  on  the  part  of  the  user  but  can  also  give 
him  more  control  in  solving  the  problem  and  perhaps  more  insight  as  to  what 
is  happening.  lf\  solving  the  associated  unconstrained  problems,  the  user  has 
the  option  of  selecting  a  direct  search  method,  thereby  eliminating  the  need 
for  differentiable  or,  in  some  problems,  continuous  functions.  Special-purpose 
MP  routines  can  be  closely  tailored  to  fit  special  applications  via  these 
techniques. 


INTEGER  PROGRAMMING 

In  addition  to  the  standard  linear  programming  proolem,  there  are 
several  special  programs  under  the  heading  of  integer  linear  programming 
(ILP).  The  problem  statement  is  similar  except  for  constraints  placed  on 
the  variables: 

Minimize  f(x)  -  CjXj  +  02^2  +  .  . . 
subject  to 

g,(x)  =  a,jx,  -hai2X2+...  +  ajnXn<bj 

8m(*)  =  ^ml’^l  +  ^m2'‘2  •"  •  •  •  +  Wn  ^  (8) 

Xj  >  0  and  xj  is  an  integer. 

There  are  further  distinctions  within  ILP  -  pure  integer,  mixed  integer,  and 

|0,l|.  The  above  problem  is  a  pure  integer  problem;  if  X] , . .  . ,  X|^  are  required 

to  be  integral  and  Xjj^.j, .  . . ,  Xj^  are  not  necessarily  integral,  then  we  have  a 

mixed  integer  problem;  and  finally  if  Xj  can  equal  only  zero  or  one,  we  have  a 

jO,l}  ILP.  Integer  programming  is  in  its  infancy,  and  some  methods,  although 

they  theoretically  exhibit  finite  convergence,  have  not  been  computationally 

successful.  The  three  ILP  routines  available  at  NELC  are  “Zero-One  Integer 

Programming  with  Heuristics, BBMIP,^^  and  OPTALG.  Zero-one  and 

BBMIP  are  SHARE  routines  which  have  been  checked  out  on  the  360/65  but 

not  tested  extensively.  BBMIP  (a  mixed  integer  routine)  has  been  tried  on  a 
in 

series  of  test  problems  and  compared  with  other  routines;  however,  the 
other  codes  are  machine-dependent,  so  the  comparison  is  not  meaningful. 

OPTALG  is  a  bound  and  scan  pure  integer  programming  code  which 
has  solved  large  problems  successfully.  It  was  developed  at  Stanford  by 
F.  S.  Hillier,  and  is  currently  operational  at  NELC.  The  routine  requires  336k 
of  core  and  a  solution  to  the  associated  linear  programming  problem;  that  is, 
we  simply  drop  the  restriction  that  Xj  be  an  integer.  This  solution,  together 


with  tlie  LP  basis  matrix  and  a  guess  at  an  initial  feasible  solution  are  then  used 
as  data  for  OPTALG.  If  the  problem  is  large  (a  maximum  of  61  rows  and  61 
columns),  then  this  data  preparation  can  be  tedious  if  done  by  hand.  The 
procedure  has  been  somewhat  automated;  the  exact  details  are  in  the  user’s 
guide  (Appendix  2). 

The  user  should  be  cautious  when  attempting  to  solve  an  ILP  with  an 
integer  programming  code,  since  it  is  possible  for  a  routine  to  solve  some  prob¬ 
lems  and  not  others,  even  if  they  are  the  same  dimension  and  fairly  similar. 
Matching  the  routine  to  the  problem  is  still  an  art.  Progress  is  being  made  in 
this  area,  but  it  will  be  some  time  before  methods  are  available  to  solve  a  general 
ILP.  A  selection  of  engineering  applications  is  presented  in  Appendix  1 . 

The  following  graph  (fig.  4)^^  gives  an  idea  of  the  sizes  of  mathematical 
programming  problems  that  can  currently  be  solved.  Tlie  abscissa  represents 
the  sum  of  both  the  number  of  constraints  and  the  number  of  variables;  integer 
programming  methods  are  still  too  inconsistent  to  include. 


Figure  4,  Magnitudes  of  currently  solvable  MP  problems. 


DESIGNING  WITH  MATHEMATICAL  PROGRAMMING  AS  AN  AID 


The  engineer  can  effect  a  better  design  in  some  areas  in  less  time  and 
at  lower  cost  with  mathematical  programming  techniques  than  with  classical 
methods.  The  key  to  success  is  in  the  word  “aid.”  for,  in  using  the  methods 
effectively,  the  engineer  must  have  command  of  his  own  field  and  understand 
some  basic  principles  of  MP  techniques.  It  is  recognized  that  typically  the 
user  of  MP  techniques  is  concerned  with  obtaining  results  quickly  without 
lengthy  excursions  into  numerical  analysis  or  mathematical  programming. 
Commercial  routines  -  ECAP,^^  MATCH, etc.  —  are  geared  to  such  a  user; 
however,  blindly  accepting  the  output  of  such  routines  can  be  disaotrous.^^ 

Also,  these  off-the-shelf  routines  are  procrustean  -  they  generally  do  not  lend 
themselves  to  modification  for  a  special  case.  General  mathematical  programming 
computer  codes  can  be  modified  for  particular  requirements.  For  example,  if 
filter  design  by  MP  techniques  is  a  frequent  task,  then  the  computer  code  can 
be  modified  to  incorporate  the  automatic  scaling  of  variables. 

In  these  remaining  sections,  we  discuss  how  an  engineer  would  proceed 
from  start  to  finish  in  using  MP  as  a  design  aid.  We  also  discuss  some  special 
topics,  such  as  duality,  that  do  not  apply  in  all  instances,  but,  if  used  properly, 
can  save  time  and  money,  both  in  setup  analysis  and  in  obtaining  a  numerical 
solution. 


FORMULATING  A  MATHEMATICAL  PROGRAMMING  PROBLEM 

In  many  design  problem  statements  there  is,  instead  of  a  single  goal,  a 
collection  of  specifications  to  be  satisfied.  This  gives  the  engineer  several 
degrees  of  freedom  in  formulating  an  associated  MP  problem.  He  has  the  option 
of  merely  satisfying  all  the  design  specifications  tthis  is  equivalent  to  finding 
an  initial  feasible  point)  or  of  singling  out  one  distinguished  requirement  and 
using  it  as  the  objective  function.  For  example,  in  reference  1 ,  a  tunable 
bandpass  filter  was  to  be  designed  which  would  satisfy  the  following  (among 
other)  specifications: 

1 .  At  the  tuned  frequency  F,  the  “insertion  loss”  of  the  filter  was  to 
be  less  than  2  dB. 

2.  At  10%  on  either  side  of  F,  the  “roll  off’  was  to  be  at  least  40  dB. 

Either  of  the  above  options  for  formulating  an  MP  problem  would  have  been 
satisfactory.  The  latter  method  was  chosen.  The  insertion  loss  of  the  filter 
was  selected  to  be  minimized  and  the  roll  off  requirements  were  treated  as 
constraints. 

To  properly  formulate  an  MP  problem,  the  user  must  explicitly 
determine  the  following; 


1 .  The  objective  (performance,  tolerance,  etc.)  that  is  to  be  accom¬ 
plished. 

2.  The  mathematic  al  relations  that  govern  the  interaction  between  the 
independent  design  variables. 

3.  The  bonnds  and  limitations  on  the  values  of  the  components  that 
guarantee  a  realizable  design. 

With  this  information  in  hand  the  designer  can  select  which  MP  approach  to 
use.  However,  since  care  is  required  in  choosing  the  objective  function  and 
providing  a  code  for  its  numerical  evaluation,  we  present  some  general 
examples  and  guidelines. 

The  objective  function  can  be  defined  as  a  measure  of  the  merit  or  the 
desirability  of  a  solution  to  a  problem,  and  its  magnitude  typically  represents 
cost,  profit,  performance,  quality,  etc.,  or  a  combination  cf  these.  The  case  of 
the  single  well  defined  objective  function  generally  poses  no  problems;  it  is 
the  combination  of  goals  which  can  lead  to  difficulties.  The  following 
examples  illustrate  various  treatments  of  multiple  goals. 

Suppose  that  it  is  desired  to  minimize  both  the  insertion  loss  of  a  net¬ 
work  denoted  by  fj(x)  and  the  cost  of  the  components  denoted  by  f2(x). 

Then  one  formulation  of  an  objective  function  f  would  be 

Minimize  f(x)  =  fj(x) -t  r*f2(x)  (9) 

where  c  is  an  appropriately  chosen  scaling  factor. 

Another  possible  choice  of  f  would  be 

Minimize  fix)  = -f2(x)/fj(x)  (10) 

Note  that  no  scaling  factor  is  required  and  the  dimensions  of  the  objective 
function  are 

-  cost($)/(unit  power  loss)  =  cost($V(unit  power  gain)  (11) 

The  next  example  illustrates  treating  secondary  design  goals  as 
constraints.  Suppose  that  high  reliability  (f3(x))  is  desired  but  the  primary 
goal  is  a  design  of  minimum  costs(f2(x)).  A  constrained  formulation  would  be 

Minimize  f2(x)  (12) 

subject  to 

-(^(x)  +  a  <  0 

where  a  is  a  tolerance  on  reliability  (mean  time  to  failure). 

A  possible  numerical  pitfall  is  combining  design  goals  in  a  haphazard 
manner.  Suppose  we  have  two  performance  indicators  Uj(x)  and  U2(x),  and 
we  desire  to  maximize  Uj  and  minimize  03  in  the  same  design.  Since  maxi¬ 
mizing  u  j  is  equivalent  to  minimizing  -u  j ,  a  possible  objective  function 
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would  be 

Minimize  f(x)  =  W2*U2(x)  -  Wj  *uj(x)  (13) 

where  the  constants  Wj  and  W2  are  required  to  make  the  function  dimensionally 
correct.  These  weights  are  important,  since  the  magnitudes  of  -u  j  and  U2  at 
the  optimum  x*  can  be  quite  different.  For  example,  if  -U|(x*)  »  0  and 
U2(x*)  *=»  1000,  then  the  effect  of  uj  is  obliterated  by  U2. 

The  above  discussion  of  objective  functions  is  intended  to  be  general. 
Aoki"^”  presents  many  detailed  examples  of  engineering  applications  together 
with  ample  background  material. 

Some  care  should  be  taken  in  the  numerical  evaluation  of  the  problem 
functions  and  their  gradients.  Coding  the  functions  offers  an  opportunity  for 
considerable  analysis  and  clever  programming.  This  task  is  done  only  once, 
against  the  many  times  that  the  functions  are  evaluated  during  the  optimiza¬ 
tion.  A  seemingly  innocent  equation  or  a  naive  way  of  combining  terms  can 
lead  to  poor  numerical  results.  The  objective  function  can  be  complicated,  and 
a  single  e/aluation  for  a  set  of  parameters  can  involve: 

1 .  A  sohition  of  a  system  of  differential  equations 

2.  Inverting  a  matrix 

3.  Table  lookups  or  interpolation 

4.  All  of  the  above 
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Gear‘d  and  Calahan"^®  discuss  similar  numerical  problems  and  some  useful 
techniques  in  applying  MP  to  engineering  design.  One  of  the  goals  of  mathe¬ 
matics  as  applied  to  computer-aided  design  is  to  free  the  applications-oriented 
user  from  the  standard  numerical  worries.  For  example,  the  routine  should  be 
able  to  analyze  the  problem  and  choose  the  best  method  for  integrating  a 
system  of  nonlinear  differential  equations.  State-of-the-art  techniques  do  not 
meet  this  goal;  thus,  it  is  still  up  to  the  user  to  make  the  proper  selection  of 
numerical  methods. 

Duality  is  a  case  in  which  careful  analysis  can  pay  off  generously.  Full 
details  with  examples  are  presented  in  the  section  on  duality,  but  we  present 
this  idea  here  to  point  out  some  analysis  and  numerical  considerations.  In 
some  applications  of  linear  programming,  problems  occur  which  have  a  much 
larger  number  of  constraints  than  variables.  To  be  specific,  we  may  have  50 
constraints  and  10  variables.  To  obtain  a  solution  to  the  LP  problem  posed  in 
tins  way,  an  LP  code  would  essentially  invert  a  50-by-50  matrix.  This 
inversion  can  be  time-consuming  and  perhaps  inaccurate  for  such  a  large 
matrix.  We  show  in  Duality  in  MP  how  to  cast  this  problem  as  a  dual  linear 
programming  problem  which  would  require  only  a  lO-by-10  matrix  to  be 
inverted.  Duality  is  an  excellent  example  of  a  little  analysis  saving  a  great 
deal  of  time  and  effort. 


Suppose  that  the  design  problem  is  clearly  stated  as  an  MP  problem. 

In  order  to  solve  the  resulting  problem,  it  may  be  necessary  to  transform  it  into 
an  equivalent  MP  problem.  The  possible  modifications  can  occur  in  the  light 
of  the  following  questions: 

1 .  Is  the  necessary  computer  code  available? 

2.  Would  duality  aid  in  obtaining  numerical  results? 

3.  Is  an  initial  point  available? 

4.  Can  the  derivatives  be  easily  calculated? 

5.  Would  a  postoptimai  analysis  be  useful? 

6.  Can  significant  benefit  be  gained  by  scaling  the  variables  or  clever 
coding  techniques? 

It  does  little  good  to  cast  a  design  problem  as  an  MP  problem,  if  we  lack  the 
numerical  means  to  solve  it.  Thus,  the  user  should  be  prepared  to  transform 
his  problem  into  a  solvable  form  or  into  a  more  useful  f  irm.  In  the  ensuing 
pages  we  address  ourselves  to  problem  modifications  «  hich  we  found  useful 
in  practice. 

SUMT  AND  CONSTRAINT  TRANSFORMATION 

The  most  common  mismatch  is  that  of  a  constrained  problem  to  be 
solved  and  a  routine  for  solving  only  unconstrained  problems.  The  sequential 
unconstrained  minimization  techniques  (SUMT)  developed  by  Fiacco  and 
McCormick^^  transform  a  constrained  MP  problem  into  an  equivalent  sequence 
(in  terms  of  having  the  same  solution)  of  unconstrained  problems.  The 
transformation  takes  place  with  the  aid  of  an  unconstrained  auxiliary  function 
which  has  the  following  form: 

m 

0(x,  r|j)  5  f(x)  +  p(r,j)  ^  G(gj(x))  ( 1 4) 

i=1 

where  rj^  is  a  parameter,  G(y)  is  a  monotonic  function  of  y  that  behaves  in 
some  well  chosen  manner  at  y  =  0,  and  p(r)  is  a  function  of  r  which  depends 
on  the  choice  of  G.  Typical  choices  require  G(y)  >  0  for  y  >  0  and  G(y)  =  0 
for  y  <  0,  or  require  that  G(y)  approach  «>  as  y  approaches  0  through  values 
less  than  zero.  The  first  choice  of  G  is  usually  associated  with  procedures 
that  are  not  concerned  with  constraint  satisfaction  except  at  the  solution 
(exterior  methods);  and  the  second  choice  of  G  is  associated  with  procedures 
which  enforce  constraint  satisfaction  throughout  the  minimization  (interior 
methods).  The  basic  idea  of  SUMT  is  the  following.  Let  x*  be  a  solution  to 
problem  I ,  that  is,  we  assume  x*  minimizes  f(x),  subject  to  gj(x)  <  0  for  i=  I , 

,  .  m.  Then  under  appropriate  coi  itions^*^  on  the  problem  functions  the 
following  theorem  holds. 


If  jxj^  j  1  is  a  sequence  of  points,  each  of  which  minimizes 
0(x,  fjj)  where  rj^  is  a  sequence  of  points  tending  to  zero,  then  we  have  the 
limit  Xi,  =  x*,  or  in  some  cases  limit  f(xi.)  -  ffx*).  So  a  constrained  problem 

\i-*oo  k-Hx. 

is  replaced  by  an  equivalent  sequence  (in  the  sense  of  having  a  common 
solution)  of  unconstrained  problems.  Common  choices  for  G(y)  are  -1/y, 

(min  [y,  0]  e  >  0,  and-log  (-y)  with  p(r)  =  r,  1/r,  r  respectively.  For  the 

first  form  of  G(y),  <^(x,  rj^)  becomes 

m 

<?'(x,r|^)=  f(x)  +  r,j^|-l/gj(x)|  (15) 

*  *  then  if 

we  have  a  feasible  point  x°  and  seek  to  minimize  0(x,  rj^),  the  term  S  l/-gj(x) 
keeps  intermediate  test  points  in  the  constraint  set.  For,  if  the  routine 
attempts  to  leave  the  feasible  set,  it  must  cross  the  boundary;  this  causes 
2  t/-gj(x)  to  approach  +  «>,  and,  since  we  are  minimizing,  the  routine  auto¬ 
matically  avoids  points  which  yield  large  values  of  0(x,  rj^).  As  rj^  approaches 
zero,  xjj  is  approaching  x*,  and  since  xj^  is  used  as  the  initial  point  to  find 
Xjj+j,  each  successive  minimization  requires  fewer  iterations.  In  Lagrange 
Multipliers  as  a  Design  Aid  we  give  some  examples. 

Variations  of  the  above  techniques  can  be  applied  in  many  situations 

where  there  is  a  mismatch  between  computer  code  and  mathematical 

O') 

programming  probiem.  For  example,  Rosen’s"^ gradient  projection  method 
solves  the  following  problem; 

Minimize  fix) 

Subject  to 
n 

S  aijXj<bi 

j=l 

Now  suppose  the  MP  problem  we  have  on  hand  has  some  nonlinear  constraints 

gj(x)  <  0,  i  =  1 . 8,  as  well  as  some  linear  ones  We  modify  the  problem 

as  follows: 

C 

Minimize  fix)  -  rj^  ^  l/gjlx)  (17) 

i=l 

n^ 

subject  to  ^  ajj  Xj  <  bj  i  =  1, .  .  .,  m 

j=l 

Thus,  we  take  advantage  of  a  good  routine  which  our  problem  does  not 
quite  fit. 


(16) 


i  =  1 , 2, . . .,  m 
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Another  method  which  allows  us  to  make  use  of  an  unconstrained 
computer  code  is  a  straight  constant  transformation.  Many  times  in  Problem 
1  the  constraints  are  fairly  simple  -  for  example,  aj  <  Xj  <  bj  or  Cj  <  a|Xj 
+  bj  <  dj  -  and  the  objective  funct'on  is  very  complicated  with  a  difficult 
gradient  to  compute.  The  simple  nature  of  the  above  constraints  adows  in 
initial  feasible  point  to  be  obtained  immediately,  and  the  constrainv  trans¬ 
formation  keeps  the  intermediate  points  feasible.  This  transformation  (as 
well  as  SUMT)  allows  a  direct  search  method  to  be  used  for  optimizing 
complicated  objective  functions.  However,  as  the  examples  show,  there  is  a 
point  of  diminishing  return  in  trying  to  transform  all  the  constraints,  since, 
as  the  constraints  become  complicated,  it  can  be  impossibie  to  find  a  well 
mannered  constraint  transformation  to  do  the  job. 

If  we  require  a  <  x  <  b,  then  the  following  transformations  keep  x 
within  this  range: 

b+a  (b-a) 

1 .  X  =  —  + - sm  y 

2  2 

2.  X  =  a  +  (b-a)  sin'^y 

eV 

3.  X  -  a  +  (b-a) -  for  y  unconstrained. 

eV-He-y 

For  one-sided  boundaries  -  that  is,a  <  x  -  we  have: 

4.  X  =  a  +  ey 

5.  X  =  a  +  lyl 

6.  X  =  a  +  y^ 

The  next  types  of  constraints  are  linear  inequalities  -  for  example, 

d  <  ax  +  b  <  c.  However,  this  is  equivalent  to  a  boundary  inequality;  that  is, 

d-b  c-b 
—  <  X  <  —  for  a  >  0. 
a  a 

If  we  have  two  linear  inequalities  in  two  unknowns,  a  transformation 
is  still  possible: 

a  <  b|  jX[  +  b|2X'j  <  c 

e  <b2ixi  +  b22X2<f 

then  let  y i  “  h|  |Xj  +  b|2X2 

y2  =  b2ixi  +b22X2 

then  X|  =  (b22yi -b|2y2)/D 

X2  =  (b2iyi  -  b]  iy2)/f> 
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where  ^  ~  *^22^1 1  "  *^12^21 
then  let  yj  =  a  +  (c-a)  sin  ^  zj 
y2  =  e  +  (f-e)  sin  ^ 

Thus,  we  have  made  Xj,  X2  functions  of  the  unrestricted  variables  Z|,  Z2,  and 
they  still  satisfy  the  inequality  constraints. 

OTHER  EXAMPLES 

1 .  Suppose  we  require  0  x  j  <  X2 

then  consider  xj=yj 

-  2^  2 
''2  =  yi'^y2 

2  2  2 
X3  =  yi+y2  +  y3 

2.  This  example  was  presented  in  Box^^ 
maximize  f  =  [9  -  (x  j  -  3)^] 


subject  to 

0<xj 

0  <  X2  Xj/v^ 

0 < X|  +\/3  (X2) ^ 6 

Initial  point 

xi  =  1,  X2  =  .5 

f=  .01336 

optimum 

xi  =  3,  X2  =Vf 

f=  1.0 

TRANSFORMATION  OF  CONSTRAINTS 
yj  =  xj  +VTx2 
y2  =  x2-xi/>/r 

Then  if  X2  -  xj/  y/3^  ^  we  have  0<Xj  +  X|  <6orxj  <3.  This  then  gives 
bounds  for  the  variables  yj  and  y2. 

-Xj/VT  <X2-Xj/\/3  <0 

=-3ls/J  <  y2  <  0.  Thus  we  have 
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0  <  yj  <6 

-y/J  <y2<0. 

This  implies 

xj  =  1/2  (yj  -\/3y2) 

X2  =  1/2  (y|/vT  +  y2) 

Thus,  we  have  the  unconstrained  problem 
maximize  f(x j(zj,  Z2),  X2(Zj,  Z2)) 
where  Xj  =  1/2  (yj  -  \/3y2) 

X2=  1/2  (yj/\/3  +y2) 

yj  =  6  sin  "^Zj 

y2  =  -  VT  +  sin  ^Z2 

The  necessary  FORTRAN  code  to  evaluate  this  objective  function,  if  we  were 
using  NFLC  DIRECT  to  solve  the  minimization  problem,  is  as  follows: 

SUBROUTINE  FN(Z,  N,  F) 

DIMENSION  Z(2) 

Y1  =6.0*(SIN(Z{1)))^  Z 

Y2  =  -2.0*SQRT(3)*(HSIN(Z(2)))**2) 

XI  =  .5*(Y1-SQRT(3)*Y2 

X2  =  .5*(Y1/SQRT(3)+Y2) 

F  =  (9-(Xl-3.)**2)’^(X2**2)/(27.0*SQRT(3)) 

RETURN 

END 

Tlio  best  transformation  to  use  depends  both  on  the  problem  and  the 
minimization  scheme  used.  For  example,  transformation  5  is  not  differen¬ 
tiable  at  y  =  0;  hence,  a  gradient  method  would  not  be  valid  at  this  point. 
Method  3  would  make  coding  the  gradient  of  the  objective  function  lengthy 
and  time-consuming  to  debug.  These  problems  do  not  occur  if  a  direct 
search  method  is  used. 


GRADIENT  APPROXIMATION 


Let  us  now  suppose  the  problem  is  matched  to  the  routine  and  further 
suppose  the  computer  code  is  a  gradient  method.  The  user  must  supply  a 
subroutine  or  function  that  will  compute  the  gradient;  that  is,  he  must 
provide  !he  following  vector: 


VfT(x)  = 


^  ^  iL\ 

dx,’dx2’  “’dx^j 


(18) 


Often,  the  complexity  of  the  objective  function  f(x)  or  of  some  of  the 
constraint  function  is  such  that  the  gradient  cannot  be  computed,  or  the 
time  required  to  derive  Vf  analytically  is  excessive.  Then,  the  designer,  in 
order  to  use  the  gradient  methods,  can  decide  to  approximate  the  gradient. 

In  the  process,  the  function  must  be  evaluated  several  times  in  the  vicinity  of 
a  point  x  for  each  approximation.  Thus,  a  trade-off  situation  arises,  since  th  e 
more  accurate  estimate  will  require  many  function  evaluations  which  may 
ultimately  be  more  costly  (in  computer  time)  than  direct  calculation  of  the 
gradient  (in  manhours).  The  problem  becomes  one  of  selecting  the  most 
accurate  approximation  for  the  least  number  of  function  evaluations. 

We  list  some  of  the  standard  schemes  for  approximating  the  derivative 
of  a  function  together  with  the  corresponding  error  estimates  (Hildebrand^®). 
Let 


fj  =  f(xi,  X2, . .  .,Xj  + A,  ...,  x„) 

(19) 

fo  =  ffxj,  X2,  . . .,  xj,  . .  .,x„) 

(20) 

f-l  =f(Xj,  ...,  Xj-A,  ...,Xj,) 

(21) 

3^f 

e  = - (xj,  ...,  Xj+t  . .  .,x„)where  |^1<A 

3x^ 

'  eA2 

3f , /3xj  =  ( 1 /2A)(-3f.i  +  4fo  -  f  1 )  + -— - 

(22) 

(23) 

2 

3fo/3xi  =  (l/2A)(f., +  f,)-^- 

O 

(24) 

3f.,/3xi  =  (l/2A)(f.,  -4fQ+3f,)  +  — 

(25) 

We  note  that  the  best  estimate  for  df/3xj  is  equation  (29),  for  if  |d^t/3x?i 
<  M3  as  ^  varies  over  |  ^  |  <  A  we  have  the  maximum  error  |  e  max  |  =  M3/6. 
It  is  interestir  that  the  point  in  question  does  not  appear  in  the  formula, 
even  though  fQ  is  generally  available.  This  additional  information  does  not 
improve  the  estimate.  The  next  question  is,  how  should  A  be  chosen? 
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Hildebrand  gives  an  optimal  A  which  is  derived  analytically  via 
3-’r(x)/9x. .  Unfortunately,  this  expression  is  generally  not  available  when  we 
are  attempting  to  approximate  3f(x)/3xj.  Also,  this  optimal  A  is  a  function 
of  the  test  point  x,  which  changes  many  times  in  the  course  of  a  minimization. 
If  this  approach  is  taken,  then  an  educated  guess  will  have  to  be  made  for  A. 
Experiments  have  indicated  that  for  the  most  efficient  operation,  A  should  be 
changed  automatically  by  the  code.  A  method  along  this  line  is  Stewart ’s'^^ 
modification  of  Davidon’s^^  minimization  routine. 

Steward’s  method  is  dependent  on  the  information  generated  by  the 
minimization  part  of  the  routine.  Using  this  information,  he  is  able  to  select 
a  good  Aj  and  initially  estimate  the  ith  component  of  the  gradient  by  the 
formula  (f  j  -  IqI/Aj.  When  this  simple  scheme  begins  to  fail,  the  routine 
automatically  switclies  to  the  central  difference  method  for  a  more  accurate 
estimate.  Aj  is  also  suitably  modified.  His  method  required  1 63  function 
evaluations  to  find  the  minimum  of  the  Rosenbrock  function  to  within  five 
decimal  places.  This  compares  with  325  for  the  direct  search  method 
(ZANGWL)  and  71  for  the  gradient  method  (CNJGAT),  to  minimize  to 
approximately  the  same  accuracy.  Unfortunately,  Steward’s  method  was 
coded  in  a  language  for  the  CDC  1604  and  has  not  been  modified  for  the 
NELC  IBM  360. 

The  final  technique  for  general  numerical  differentiation  which  we 
discuss,  is  a  code  from  the  IBM  Scientific  Subroutine  Package^  called 
DDCAR.  DDCAR  uses  an  extrapolation  technique  to  obtain  highly  accurate 
estimates  of  the  gradient.  This  method  does  not  need  an  optimal  A  to  give 
good  results.  Howevt-,  we  must  pay  the  price  of  using  a  larger  number  of 
function  evaluations  than  the  central  difference  method  would  requirs,  for 
each  estimate  of  Vf.  A  code  and  results  of  using  DDCAR  to  compute  the 
gradient  in  minimizing  the  Rosenbrock  function  with  the  gradient  method 
NELC  CNJGAT  are  presented  in  table  1.  1850  function  evaluations  were 
required  to  obtain  these  results. 

The  methods  presented  above  are  general  and  can  be  used  for  a  wide 
class  of  functions.  Some  applications  may  lend  themselves  to  special  methods 
for  estimating  the  gradient.  For  example,  Calahan"^"  discusses  methods  for 
numerically  evaluating  the  gradient  in  a  network  optimization  scheme.  The 
method  is  highly  specialized  for  this  type  of  problem.  He  uses  a  “calculus  of 
variations’’  approach  together  with  numerical  integration.  The  section  on 
gradient  calculation  is  a  good  example  of  the  success  of  careful  analysis. 


IMPLICIT  RE/VL*a  (A-H,0-ZI 
COMMON  KOUNT 

DIMENSION  X(3C»»  H(3C,30»,  G(30J,  S1301,  S IGMA ( 30 » ,XX ( 20 ) 

EXTERNAL  FG  BOX 

X(  1)--1.2D00 

X(  2)--1.0D00 

KOUNT=  0 

NMAX  =  30 

MPPNT  =  1 

N=  2 

I  START  =  N  ♦  1 
FPSLON  =  1. 00-10 
FEST  =  0.0 
ITEPBD  =  80 
SBOUND  =  1.0D-C9 
ICNT  =  0 
IRE  ST  =  0 

iter  =  0 

NMl  =  N  -  1 

DO  30  I  =  1,  NMl 

H(I  ,I  )  =  KO 
IPl  =  I  +  I 

DO  30  J  =  IPl,  N 

H( J, n  ^  0.0 
30  H( I.J  I  =  C.O 

H(N,N)  =  1.0 
DO  2C  J=l,20 
20  S(J)  =  0. 

CALL  CNJGAT(FGBOX,N»NMAX,X  ♦  I  TERBD , EPSLON, F EST , MPRNT  , 
ISTART,ITE  1  P,F,G,SI 
PRINT  qqq^KOUNT 
9QO  FORMAT!  1 X , • KOUNT • , I  5 ) 

END 

FUNCTION  USERF  (X,N) 

IMPLICIT  RFAL*8  (A-H,0-Z) 

DIMENSION  XIN) 

T1  =  X(2)-X(  1)*X(  1  )*X(  1) 

T2=X( 1) *1.0000 

USERF=100.CO0C*Tl*Tl  +  T2*T2 

RETURN 

END 

FUNCTION  FUNC(H) 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  KOUNT/r,P.An/XX,IVAR,NN 
DIMENSION  XTEMP(20)  ,XX{20) 

DO  30  J  =  1,NN 

IF  ( J  .EQ.  IVAR  J  GO  TO  29 

XTEMP(J)  =  XX(J> 
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Gfl  Tn  3G 

XTEMP  (IVAR)  =  XX(!VAR)  +  H 

CONT INUE 

FUNC=USF‘’F  (  XTt 

KCIUMT  =  KOUNT  +  1 

PFTUPN 

EK'D 

SUPPPHTINL  FG  PDXC  X  fF  t  Gt  N) 
IMPLICIT  RFAL*P  (A-H,n-Z) 
EXTERNAL  FUNC 

DIMENSION  XIDfGd)  »XX(2C) 
COMMnN/GRAO/XX, I VARfMN 
NN  ^  N 

DO  10  J=lfN 
XX (J  »  =  x{  J ) 

IVAP  =  1 

F  =  FUNCIOI 

DO  20  K-lf  N 
IVAR  =  K 

WZ  =  .001 

CALL  DUCARIO. fWZ  tl»FUNC»Y) 

G{K I  =  Y 

RETURN 

END 
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DUALITY  IN  MP 


In  linear  programming  and  certain  SUMT  methods  'v.'th  convex 
functions,  the  possibility  of  using  duality  to  gain  computational  efficiency 
should  be  considered. 

Duality  occurs  in  many  areas  —  mathematics,  engineering,  economics 
physics,  etc.  In  a  mathematical  or  engineering  context  it  implies  that  two 
concepts  or  systems  have  a  specific  mathematical  relationship.  We  give  an 
example  from  circuit  theory  before  proceeding  with  duality  and  MP. 

Consider  the  following  series  RCL  network  (fig.  5).^^ 


Figure  5.  Dual  circuits. 


The  mathematical  equation  representing  (A)  is 
di  1 

L  —  +  Ri  +  —  J  idT  =  e(t) 
dt  C  ^ 

If  we  interchange  e  and  i  and  L  and  C  and  replace  R  by  G  =  1/R,  we  have 

_  de  ’  r  ^ 

C  —  +  Ge  +  —  /  edr  =  i(t) 

dt  L  -^0 

which  is  the  equation  which  models  the  parallel  network  (B).  Thus,  an 
equation  of  the  form 


(26) 


(27) 


dx 

—  +  bx  +  c  /  xdT  =  y(t) 
At  Jq 


dt 


can  represent  either  circuit  (A)  or  circuit  (B)  depending  on  the  parameters 
a,  b,  that  is,  it  has  a  dual  function,  and  we  say  that  (A)  is  the  dual  of  (B). 

In  linear  programming,  we  have  the  most  straightforward  application 
and  complete  theory  of  duality.  Let  x  and  c  be  vectors  of  length  n,  b  a 
vector  of  length  m,  and  A  an  m-by-n  matrix.  Then  the  linear  programming 
problem  seeks  to  find  the  vector  x  which  produces 


a  minimum 

T 

z  =  C  X 

subject  to 

Ax=^b 

(28) 

x>0. 

To  this  linear  programming  problem  there  corresponds  an  associated  problem. 
Maximize  v  =  b^w 

subject  to  (29) 

T 

A^w  >  c 

w>0 

where  w  is  an  m  vector.  Problem  (28)  is  called  the  primal  problem  while 
problem  (29)  is  referred  to  as  the  dual  problem. 

In  the  primal  formulation  we  have  m  constraints  and  n  variables,  and 
just  the  opposite  in  the  dual  formulatic;  Some  facts  about  these  corre¬ 
sponding  problems  which  are  pertinent  to  this  section  are; 

1 .  The  primal  problem  has  a  bounded  solution  x*  if  and  only  if  the 
dual  has  a  bounded  solution  w*. 

2.  w*  can  be  obtained  from  x*  and  vice  versa. 

3.  The  dual  of  the  dual  is  the  primal. 

4.  In  the  case  of  bounded  solutions  we  have  z*  =  v*. 

The  efficiency  of  linear  programming  computer  codes  decreases  as  the 
number  of  constraints  increases.  It  is  this  trait  which  makes  the  study  of  the 
duality  relationship  worthwhile  from  a  computational  standpoint.  The 
following  example  shows  how  this  loss  of  efficiency  can  be  lessened  by 
judicious  use  of  the  dual. 

Consider  the  following  programming  problem 

Minimize  z  =  -4x  j  -  3x2 

subject  to 

Xj  <  6 
X'^  ^  8 

Xj  +  Xt  <  7  (30) 

3xj  +  X2  <  15 
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-X2<1 

Xj  >0, 

Note  that  X2  can  be  negative.  To  pose  this  as  a  standard  ^LP  problem,  we 
replace  the  variable  X2  by  the  difference  of  two  positive  variables;  that  is,  we 
let  X2  =  X2  -  X2 .  Now  the  problem  can  be  written  in  standard  form. 

Minimize  -4x  j  -  Sfx^  -  x^’) 
subject  to 

xj  <6 
X2-x^'<8 
X  j  +  X2  ■  ^2  ^  7 

3xj+x^-x^'<15  (31) 

-x^  +  x^'<l 

X|,  xi),  X2  >  0 


Writing  the  constraints  in  matrix  notation,  we  have 


■*  * 

1  0  0 

6 

0  1  -1 

’^l 

8 

1  1  -1 

^2 

< 

7 

3  1  -1 

^2 

15 

[0  -1  1 

W  a 

1 

In  solving  this  LP  problem,  we  add  another  positive  variable  xsj,  called  a  slack 
variable,  to  each  row,  to  make  each  inequality  an  equality.  The  problem  then 
becomes: 


^MPS/360  accepts  unrestricted  variables  and  automatically  makes  this  transformation. 
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Minimize  -4xj  -  3x2  +  3x2  +  Oxsj  +  . .  .  +  Oxsj 
subject  to 


1  0  0  10  0  0  0 
0  1-10  10  0  0 

1  1-10  0  10  0 

3  1-10  0  0  10 

0  -1  10  0  0  0  1 


(32) 


Note  that  each  slack  variable  appears  in  the  cost  row  with  zero  for  a 
coefficient. 

The  simplex  method  of  LP  performs  an  iterative  process  on  five 
selected  columns  of  this  augmented  matrix  to  obtain  a  basis  inverse  matrix. 
This  matrix  is  then  used  to  obtain  the  solution.  It  is  this  inversion  process 
which  causes  the  efficiency  of  LP  codes  to  decrease  as  m  increases.  We  now 
cast  problem  (3 1 )  in  its  dual  formulation,  which  will  reduce  the  number  of 
rows. 

Maximize  v  =  6w  j  +  8w2  +  +  1  Sw^  +  wj 

subject  to 

f'^ll  _  . 

1  0  1  3  o']  W2  -4 

a'^w=  0  111-1  W3  >  -3  (33) 

0  -1  -1  -1  ij  W4  L^. 

'^5 

>•  « 


Now  from  each  row  we  subtract  a  positive  variable  wsj  and  write  the 
constraint  matrix  as: 
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We  now  have  a  formulation  similar  to  problem  (33);  however,  to 
solve  this  problem,  it  is  necessary  to  “invert”  only  a  3-by-3  matrix.  Even 
though  the  number  of  variables  is  the  same,  the  smaller  number  of  constraints 
makes  this  formulation  more  efficient.  If  we  were  to  solve  the  dual  formu¬ 
lation  using  MPS/360,  since  the  dual  of  the  dual  is  the  primal,  the  solution 
to  problem  (31)  would  appear  under  the  DUAL  ACTIVITY  heading. 

There  is  another  form  of  the  primal-dual  relationship  which  has  both 
formulational  and  computational  advantages;  it  is  the  unsymmetric  form. 

The  primal  problem  can  also  be  stated  as:  find  a  column  vector  x* 

which 

...  T 

minimizes  z  =  c^x 
subject  to 

Ax  =  b  (34) 

x>0 

The  original  LP  problem  (28)  can  be  put  in  this  form  by  adding  slack  variables 
to  transform  the  inequality  constraints  into  equality  constraints.  Then  the 
unsymmetric  dual  to  (34)  is 

Maximize  v  =  b^w 

subject  to  (35) 

A^w  <c 

We  notice  in  (35)  that  there  is  no  restriction  on  the  sign  of  w.  This  is  most 
useful  in  using  linear  programming  as  an  analysis  tool,  and  to  condense  the 
problem  size.  The  three  previous  properties  for  the  symmetric  form  of  the 
dual  also  hold  for  the  unsymmetric  form. 

We  return  to  problem  (30)  to  give  a  simple  example  of  the  usefulness 
of  this  unsymmetric  form.  The  problem  is  written  as 

Maximize  v  =  4w|  +  3w2 

subject  to 
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Note  that  W2  is  unrestricted  in  sign  and  that  the  last  row  of  keeps  Wj 
nonnegative.  This  is  the  unsymmetric  dual  formu’ation  of  problem  (30). 
Since  W2  is  unrestricted  in  sign,  to  solve  this  problem  via  a  standard  LP  code, 
we  should  have  to  replace  W2  by  w^  -  w^',  *is  before.  But  if  we  consider  the 
associated  primal  problem,  this  requirement  disappears;  i.e.. 

Minimize  6xj  +  8x2  +  7x3  +  ISx^  +  Xj 
subject  to 


To  solve  this  primal  problem  we  must  invert  only  a  2-by-2  matrix 
rather  than  a  5-by-5  as  in  the  dual  case.  When  we  solve  this  program  using 
MPS/360,  the  optimal  w  will  appear  as  a  DUAL  ACTIVITY. 

This  method  accomplishes  two  things: 

1 .  There  is  no  increase  in  the  number  of  variables. 

2.  If  m  is  much  larger  than  n,  then  the  primal  has  fewer  constraints 
and  will  generally  be  faster  to  solve. 

Of  course,  the  form  we  decide  to  use  will  depend  on  the  relative  sizes 
of  m  and  n. 

We  present  another  example^^  which  will  illustrate  both  the  use  of 
the  unsymmetric  dual  and  the  utility  of  LP  in  the  area  of  applied  mathematics; 
viz.,  linear  boundary  value  problems.  Consider  the  following  problem.  Find 
a  solution,  y,  to  L[y]  =  r(x)  over  [a,  b] ,  where 

n 

L(yl=^  fj(x)y^)  =  yx)y +  f|(x)y('U  . . +  fj^(x)y^")  (38) 

j=0 

with  boundary  conditions  Vj[yj  =  Pj,  j  =  1, 2, . . .,  n 

n-1  (39) 

where  Vjly)  =  ^  (oj  j^y((^)(a)  +  j^y(*^\b)). 

k-0 

and  y(‘^  is  the  ith  derivative  of  y(x)  with  respect  to  x. 
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Some  boundary  value  problems  do  not  have  a  clos''.d  form  solution  or  even  a 
solution  in  the  limit;  however,  we  should  still  like  some  information  and  in 
most  instances  an  approximate  solution  is  sufficient.  The  approximation  is 
made  in  the  min-max  sense;  i.e.,  if  f*  is  the  theoretical  solution  to,  say,  a 
differential  equation  over  an  interval  [a,  b] ,  then  an  approximation  f^^  is 
sought  to 

minimize  j  maximum  1  f^^fx)  -  f*(x)|  (40) 

x  e[a,  b] 

This  formulation  lends  itself  to  an  application  of  mathematical  programming. 
We  approximate  the  solution  by  a  sum  of  functions  over  the  interval  [a,  b] . 
To  do  this,  we  partition  [a,  b]  as  follows: 

a  <  x j  <  X3  .  . .  <  Xj^  <  b,  m  >  p  and  seek  parameters  aj , .  . .,  ap 

which  minimize 

may  I  L(y*(Xj))  -  r(xj)  |  (41) 

1  <  i  <  m 

where 

P 

y*(x)  =  yQ(x)  +  ^  ajyj(x),  (42) 

J=1 

y^ix)  satisfies  Vj.(yQ)  =  Tj.,  r  =  1, . . .,  n 
and  yj(x)  satisfy  Vj.(yj)  =  0.  j  =  1 , . . p 

Substituting  equation  (42)  into  equation  (41),  and  introducing  an  additional 
variable  e,  we  wish  to  find  aj, . . ap,  e  which  minimizes  e  subject  to 

P 

I  L  (yo(xj)  +  ^  ajyj(xj))  -  r(xj)  |  <  e 
j=l 

for  i  =  1 , 2, .  .  .  m 

For  LP  to  be  applied,  the  constraints  must  be  linear;  thus,  we  have 
Minimize  e 
(ai,  a2,  •  . 
subject  to 

P 

-e  aiL[yj(Xj)l  +  Llyglxj)]  -r(Xj)<e  j  =  1,  2. ...  m, 
i=l 


or,  rewriting  as  two  single  inequalities, 


P 

ajL(yj(Xj))- L(yQ(Xj))  +  r(Xj)<e  (43) 

P 

^  ajL(yj(Xj))  +  L(yQ(xj))  -  r(Xj)  <  e  (44) 


Collecting  all  the  parameters  on  the  left  side  of  the  inequality  yields 


- ^  ajL(yj(Xj)) Uy^fXj)) - r(Xj) 


i=l 

P 


^  ajL(yj(Xj))-t->r(Xj)-L(yQ(xj)) 
^1 


for  j  =  1, 2, ...  m 


where 


a  =  xj  <  x^  . . .  <  Xjj^  =  b  is  a  partition  of  [a,  b] . 


Form  the  p-by-m  matrix: 


G  = 


L[yi(xj)]  . . .  L[yj(Xn,)] 


Llyp(xi)]  . .  L[yp(Xn,)lJ 


the  2m  vector 

c'*'=  |L(yo(xi)-r(xj),...L(yo(Xj,^)-r(Xjj,),r(x,)-L(y^,(xi)),  .  . ., 
r(Xm)-L(yo(Xm))j 


the  (p  +  1 )  vectors 

a^=  (aj, . .  .,ap,  e) 

bT'  =  (0.  .  .0,-1) 


and  the  matrix 


Since  there  is  no  restriction  on  the  sign  of  the  parameters  aj  (i=l , . .  p),  we 

cast  the  programming  problem  in  the  dual  formulation;  i.e., 


Maximize  -e  =  «^  •  b 

(aj,a2,  ...ap, 


(45) 


subject  to 


or  more  succinctly  as 
T 

Maximize  a‘b 

subject  to  (dual) 

A^a<c 


then,  since  the  dual  of  the  dual  is  the  primal,  we  have 
T 

Minimize  c  x 

subject  to  (primal) 

Ax  =  b 

X  >  0  where  x  is  a  2m  vector. 


This  primal  problem  is  then  solved  numerically  and  the  optimal  a  vector 
appears  as  the  DUAL  ACTIVITY  vector. 

^Note  inequalities  (43)  and  (44)  keep  e  >  0. 
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In  summary,  ine  basic  advantages  are  as  follows:  We  wish  para¬ 
meters  a  j  , Up,  which  give  the  best  estimate  in  a  min-max  sense  of  the 
solution  to  the  linear  boundary  value  problem.  The  parameters  are 
unrestricted  in  sign  and  the  number  of  points  m  in  the  interval  of  solution 
is  large.  To  cast  this  problem  as  a  primal  linear  programming  problem  would 
have  required  2p  +  1  variables,  since  the  aj  are  unrestricted,  and  2m  constraints. 
By  first  writing  the  approximation  problem  as  a  dual  linear  programming 
problem  (45),  we  need  only  p  +  1  variables;  then  transforming  to  the  primal 
reduces  the  number  of  constraints  from  2m  to  p  +  1 .  We  give  two  examples 
which  use  the  above  method.  Example  1  A  homogeneous  equation 
with  inhomogeneous  boundary  conditions 

L(y)=  xy"  -(x+l)y'  - 2(x-l)y  =  0; y(0)  =  l,y(l)  =  0on  [0,  1] 

P 

y*(x)  =  yptx)  +  ^  ajyj(x) 
j=l 

with 

yQ(x)=  l-x;yj(x)  =  xj(l-x) 
j=  1,2,  ...,p 

then 

Llyg)  =  2x^  -  3x  +  3 

Uyj]  =  2x^'*’^'  +  -(j^  +  j'3)xj  +j(j-2)xj"^ 

xj  =  (i-l)(0.05).  i=l,2,...,21 

for  p  =  3,  aj  =  3.0374706 
ap  -  -0.78655970 

33  =  0.88391133 

e  =  0.03747406,  and  the  maximum  error  occurred  at  x  =  0. 

Example  2:^^  Inhomogeneous  case  with  homogeneous  boundary  conditions 
Elyl  =  y"  +(l+x^)y  =  -l,y(l)  =  y(-l)  =  0 
yj(x)  =  1  -  x“J 

L[yj]  =  ^  -  2j(2j-l )x^-i"^  +  x“  +  i 


on  [-1, 11 


Let 


z  =  x^,  Ltyj]  =  -  2j(2j  -  l)zj"^  +  z  +  1 

Zi  =  (i-1)(0.05)  i=l,2,  ..,21 

forp  =  3,  aj  =0.9675,  32“*'  a3  = -0.0285 

e  =  0.0029 

The  preceding  paragraphs  have  pointed  out  some  of  the  numerical 
advantages  of  duality.  Duality  concepts  also  have  application  in  nonlinear 
programming;  however,  so  far  these  have  been  limited  to  SUMT  methods 
involving  convex  functions.  In  this  case,  the  duality  theory  provides  a  lower 
bound  on  the  value  of  the  unconstrained  problem  to  test  for  convergence;  see 
reference  24.  Its  use  in  nonlinear  problems  is  no!  widespread. 

LAGRANGE  MULTIPLIERS  AS  A  DESIGN  AID 

Another  influence  on  choice  of  routine  or  method  to  solve  an  MP 
problem  is  postoptimal  analysis. 

When  a  constrained  optimization  problem  is  solved  via  sequential 
unconstrained  minimization  techniques  (SUMT),  additional  information  is 
available  to  the  designer  foi  a  sensitivity  analysis;  i.e.,  postoptimal  analysis. 
Suppose  X*  is  the  solution  to  the  following  MP  problem. 

Minimize  f(x) 

subject  to  (46) 

h|(x)<bj,  i=l,2, ...,  m 

The  engineer  wishes  information  as  to  how  f(x*)  will  change  if  b  is 
changed  a  little;  i.e.,  if  the  design  requirements  are  changed,  how  wi^'  the 
performance  be  affected?  The  interesting  result  is: 

=  ,47, 

ODj 

where  Xj  is  a  generalization  of  the  classic  Lagrange  multiplier  for  finding  the 
extrema  with  side  conditions.  The  proof  of  (47)  is  contained  in  reference  44. 
We  recall  from  advanced  calculus  that  an  extrema  problem  with  side 

conditions  was:  find  the  minimum  (maximum)  of  p(x  j .  X2 . Xj^)  subject 

toqj(Xj,  . .  .,  x„)=  . .  .  =  qj^(x|,  ....  x„)  =  0.  To  solve  this  problem  we 
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formed  the  Lagrangian  L(X| . ’‘n> 

then  solved  the  (n+m)  system  of  equations: 


i=l 


9L(x,  X) 
3xj 

qjtx)  =  0 


i  -  1,  . .  n 
i=  1,  .  .  m 


for  (x  I ,  .  . Xj^)  with  the  Xj’s  being  introduced  to  help  find  the  xj.  The 
Kuhn-Tucker  theorem'^^  allows  us  to  define  a  meaningful  and  useful 
Lagrangian  for  inequality  constraints. 

We  discuss  the  Kuhn-Tucker  theorem  to  allow  a  generalization  of 
Lagrange  maltiplier,  and  then  discuss  SUMT  methods  to  illustrate  obtaining 
the  Xj’s  numerically. 

Rewriting  the  constraints  to  problem  (46)  as  gj(x)  =  hj(x)  -  bj  <  0, 
we  obtain  problem  (1); 

Minimize  fix) 

subject  lu  (48) 

gj(x)  <0  i=  1,  .  .,  m. 

Tlien  the  Kuhn-Tucker  theorem  says  essentially  the  following: 

Lcl  X*  be  a  solution  to  the  above  prob’em  and  assume  the  boundary 
Oi  the  constraint  set  has  no  “cusps;”  then  the  following  conditions  hold. 


1 .  There  exist  multipliers  Xj  >  0,  i  =  I , . . m  such  that 

X,  gj(x’*')  “  0.  i=l,  2,  . .  .,  m 
n 

2.  +  ^  XjVgj(x’','  =  0 

i=l 

The  following  example  iliust.  rtes  the  Kuhn-Tu-ker  condition: 

Minimize  fi;  j ,  X2)  =  j  -  2)^  +  (X2  -  1  )^ 
subject  to 

glfxp  >  2)^  -X2  +  x|<0 
g2(x  1 ,  X2)  =  -2  +  X I  +  X2  <  0 
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The  global  solution  to  the  above  is  (x  j*,  x''‘2)  =  d .  1  )•  The  gradients 
at  the  optimum  are 


2 

1 

-2 

vgl  = 

1 

.  - 

II 

> 

1 

and  Vf  = 

0 

and  the  multipliers  are  Xj  =  X2  =  2/3.  Finally  we  can  write; 

-Vf=+2/3  vgi  +2/3  vg2- 
Geometrically  we  have  the  following  (fig.  6); 


f(i,  1) 


Vg,!!.  1) 


Figure  6.  Kahn-Tucker  conditions. 


Now  that  we  have  seen  a  need  for  knowing  the  Lagrange  multipliers 
and  have  seen  the  geometrical  interpretation,  we  turn  to  SUMT  and  the 
numerical  calculation  of  the  multipliers.  The  constrained  optimization 
problem  is  transformed  into  an  unconstrained  one  through  the  use  of  an 
auxiliary  function,  0(x).  fl(x)  does  one  of  two  things,  0(x)  -><»  as  x  0  or 
0(x)  =  0  for  X  feasible  and  6(x)  positive  for  x  -nfcasible.  The  three  forms^^ 
of  0(x)  we  discuss  are: 

m 

Barrier:  B(x)  =  l/gj(x)  (49) 

i=i 
m 

Penalty:  P(x)=  ^  (Max  (gj(x),  0))''*'^  (50) 

i=l 

m 

Logarithmic:  L(x)  = iog(-gj(x))  (51) 

i=l 

SUMT  works  as  follows:  Using  0(x)  =  B(x),  P(x)  or  L(x),  we 
transform  the  constrained  problem  into  a  sequence  of  unconstrained  problems, 
as  in  SUMT  and  Constraint  Transformation. 

Minimize  D(x,  r|^)  =  f(x)  +  p(r|^)0(x)  (52) 

for  rj^  >  0.  Then  if  x*  is  the  optimum  for  the  constrained  problem  (48)  and 
x(rjj)  ~  X|j,  the  optimum  for  problem  (52),  it  can  be  shown  that  limi*  Xj^  =x''‘; 

k-+oo 

i.e.,  (r|^-^0).  Thus,  we  replace  the  constrained  problem  by  a  sequence  of 
unconstrained  problems. 

In  the  barrier  and  logarithmic  cases  as  gj(x)  ->•  0  -  i.e.,  x  attempts  to 
leave  the  constraint  set  -  gj(x),  must  approach  0,  and  this  causes  B(x)  and 
L(x)  to  increase  rapidly.  Since  the  routine  wishes  to  minimize  D(x,  r),  and 
we  have  an  initial  feasible  point,  intermediate  x’s  are  chosen  feasible. 

In  the  penalty  function  approach  the  test  points  are  allowed  to  leave 
the  constraint  set;  however,  when  they  do  so,  a  positive  amount  is  added  to 
the  objective  function.  Again,  since  we  ,<eek  to  minimize  Dtx,  r),  points  are 
selected  within  the  constraint  set.  When  3(:U  or  L(x)  is  used,  we  have  an 
interior  point  method;  and  when  P(x)  is  used,  we  have  an  exterior  point 
method.  With  each  type  of  0(x)  a  slightly  different  technique  is  used  to 
recover  the  Xj. 
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Barrier  Method: 
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Now,  in  a  similar  argument,  let 

=  r,,/-gi(x,^)  (59) 

(k) 

Since  (-gj(x*))  >  0  in  ihe  feasible  set,  we  have  X  .  >0.  Thus,  we  can  form 

the  Lagrangian  as  before  and  equation  (59)  generates  the  multipliers. 

We  next  give  an  example  to  illustrate  this. 

Lootsma^^  solves  the  following  problem 

3 

Minimize  f(x)  =  ■  bx^"  +  I  Ixj  +  x^ 

subject  to 

g|(x)  = +x^  +  x^  -  x^  <0  f60) 

g2(x)  =  -x^-X2-X2-4<0 

g3(x)  -  +X3  -  5  <  0 

g4(x)  =  -X|  <0 
g5(x)  =  -X2  <  0 


g6(x)  =  -X3  <  0 

The  optimum  is  x*  =  (0,\/T,>/y)  with  gj(x*)  =  g2(x’'‘)  =  g4(x*)  =  0.  The 
theoretical  multipliers  are 

X,  =X2=\/278,- 0.1 767  7  66 


X4  -  1 1 ;  X3  -  X5  -  X^  -  0 
For  r|^  =  10'^,  the  numerical  multipliers  are 


(k)_ 

4 


'k 

10'^ 

-gl(Xk) 

5.66046  10'^ 

■■k 

10'^ 

-82^’‘k^ 

6.65324  10'^ 

■■k 

10"^ 

-g4(x,,) 

9.2007  10'® 

=  0.176664 

0.15030 

10.8687 


X|^  =  (9.2077  10'^,  1.4121,  1.41422) 
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The  final  example  is  the  Penalty  Function  method. 

Let 

C(x,  r)=  f(x)+-P(x), 
r 

where 

m 

P(x)  =  ^  (max.  (gj(x),  e>0. 

i=l 

Then  at  the  optimum, 

VjjC(Xk,  r^)  =  0  which  implies  Vf(x)  +  —  ’7p(x)  =  0  (61 ) 

where  VP(x)  =  +2(l  +  e)  [max.  (gj(x),  0)]^  ygj(x) 

If  we  let  max.[gi(x),  0]  ^ 

‘  'k 

(k) 

then  X  .  >0,  since  we  have  an  exterior  point  method. 

Thus,  (6 1 )  becomes  Vf(x)  +  SXj Vgj(x)  =  0, 
and,  letting  denote  the  optimum  for  rjj,  v  nave 

Xj  =  lim  X^!^^  =  lim  —  (max.  (gj(Xk),  0] 
lc-K»  *  k-»oo  ‘‘k 

In  this  way  we  can  obtain  the  Xj’s. 

In  solving  the  Lootsma  problem  by  the  penalty  method  with  e  =  1 ,  we 
obtained  the  following  results  for  r^  =  .1*10*^; 


Xk  =  (.1443*10’'^,  1.4142,  ...,1.4141,  .  .  .) 
g,(Xk)  =  .87971 10*  10-6;  g2(Xk)  =  0.8797477*  10'^ 
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F 

I 


I 

i 

\ 


X^4)=^^(max.(g4,  0))^  =  2 


j0.388414*10-l^ 

(  10-5  i 


0.7768*10"^^9fc  11 


The  last  four  constraints  of  problem  (60)  were  transformed  by  the  following 
equation  and  not  included  as  penalties: 

2 

xi  -y, 

2 

X2  =  Vs 

X3  =  5  sin^  (y(3)), 

(k) 

When  these  transformations  are  used,  X  ^  cannot  be  obtained  by  the  above 
formula. 


INITIAL  POINT  AND  SCALING 

Two  final  topics  in  solution  strategy  depend  a  great  deni  on  the 
specific  problem:  initial  point  and  scaling  of  variables.  If  the  constraints 
are  numerous  and  complicated,  then  finding  an  initial  feasible  point  can  be 
an  additional  problem  in  itself. 

Zoutcndijk'^^  gives  a  “simple”  trick  for  a  transformation  which 
replaces  the  original  problem  by  an  equivalent  problem  for  which  an  initial 
feasible  point  can  readily  be  found.  Suppose  we  are  given: 

Minimize  f(x) 
subject  to 

gj(xl<0  i=  1, 2,  . . .,  m 

and  no  initial  point  x"’  can  be  found  by  inspection  01  engineering  knowledge 
such  that  gj(x°)  <0  for  i  =  1,  2.  .  .  .,  m.  We  now  form  the  following  problem: 

Minimize  fix)  + 
subject  to  gj(x)-pj{<0 


p^>0  i  =  1,  . m 

where  p  is  a  large  number,  ^  an  additional  variable,  and  Pj  -  1  if  gjfx®)  >  0 
and  pj  =  0  if  gj(x°)  <0.  If  f  >  max.|gj(x°)  for  all  i  such  that  gj(x°)  >  o|> 
rhen  the  point  (x°,  {)  is  feasible  for  the  modified  problem.  In  reference  17 
Zoutendijk  proves  that,  forp  sufficiently  large,  the  modified  problem  will 


50 


have  the  same  solution  as  the  original.  Note  that  in  the  attempt  to  minimize 
f(x)  +  ^  is  driven  to  zero,  but  to  satisfy  the  constraints  x  must  be 

simultaneously  chosen  feasible  (for  the  original  problem).  This  is  the  part  of 
mathematical  programming  in  which  the  user’s  previous  engineering  and 
design  experience  pays  off,  since  the  better  the  initial  point,  the  faster  the 
routine  will  converge. 

Proper  scaling  of  variables  is  still  an  art.  Most  computer  routines  lor 

solving  mathematical  programming  problems  are  designed  to  follow  steep 

curved  valleys  and  sharp  ridges  in  locating  an  optimum.  For  a  one-shot 

problem,  scaling  the  variables  is  not  worth  the  time;  however,  if  a  code  is  to 

be  written  to  solve  a  large  class  of  similar  problems,  then  scaling  may  be 

worthwhile  in  long-term  savings  of  machine  time.  We  give  an  example, 
aq 

Pierre  °  which  transforms  a  poorly  shaped  objective  function  into  a  “nice” 
bowl-shaped  function  which  can  quickly  be  minimized.  Unfortunately, 
scaling  is  somewhat  limited  to  unconstrained  problems;  in  a  constrained 
problem  attempting  to  make  the  objective  function  easy  to  minimize,  the 
scaling  might  destroy  any  useful  properties  the  constraints  enjoy. 

Example: 

2  2 

If  f(xj,  X2)  =  Xj  +  10xjX2  +  lOOx^,  then  f  has  a  narrow  valley  and 
a  minimum  at  Xj  =  X2  =  0.  To  transform  f  into  a  more  desirable  shape,  we 
eliminate  the  cross  product  lOx  jX2  by  letting  Xj  =  Z]  -  abz2  and  X2  =  bz2. 
Then  f  becomes 

f=  z^  +  (-2ab  +  10b)  ZjZ2  +  (-lOab^  -t-  a^b^  +  100b^)z2  • 

Letting  a  =  5  makes  (-2ab  +  10b)  =  0  and  b  =  (1/75)*^^  solves  -50b^  +  lOOb^ 
0  2  2 

+  25b'^  =  1 ;  thus,  f  =  Zj  +  Z2,  which  is  easily  minimized  by  almost  any 
routini.^. 

It  is  highly  recommended,  when  the  variables  are  bounded  or  the 
range  of  the  variable  is  known,  that  the  vanable  be  normalized.  For  example, 
if  the  variable  Xj  lies  in  the  range. 

Li<Xi<Ui, 

then  an  appropriate  normalization  is  given  by; 

0<yi<  1.0 


where 


Normalization  is  particularly  useful  in  direct  search  algorithms  for  which  we 
must  arbitrarily  choose  a  step  size  search  increment,  A.  It  can  readily  be 
1%  of  1 .0  (0.01 ),  which  allows  each  variable  to  be  changed  proportionally, 
or  any  other  realistic  choice  depending  on  our  knowledge  of  the  problem. 

In  summary,  scaling  an  objective  function  is  generally  not  practical, 
and  we  must  rely  on  the  properties  of  the  algorithms  to  find  the  optimum; 
however,  normalization  is  recommended  wherever  possible. 


SUMMARY 

Mathematical  programming  is  a  broad  subject  with  many  varied  applica¬ 
tions.  Not  all  top'  s  in  MP  ere  applicable  to  engineering  design.  We  have  tried 
to  delineate  those  areas  wlK.'h  are  useful  and  the  corresponding  capabilities  at 
NELC.  The  computer  code;  described  have  proved  reliable  on  a  wide  range  of 
problems;  however,  MP  is  an  expanding  area,  and  the  new  algorithms  constantly 
being  developed  could  render  some  of  these  codes  obsolete.  The  applications- 
oriented  user  should  not  accept  this  list  as  complete  or  give  up  if  his  particular 
problem  does  not  match  a!;  available  routine.  Perhaps  a  literature  search  will 
yield  the  appropriate  metliod.  The  high-speed  digital  computer  has  provided 
an  impetus  to  develop  algorithms  to  solve  MP  problems  which  previously  were 
too  large  to  be  handled.  It  is  recommended  that  a  continued  effort  be  main¬ 
tained  to  keep  NELC  up  to  date  m  the  area  of  solving  MP  problems  numerically. 
Logically,  this  should  be  a  function  of  either  a  Center-wide  computer  users’ 
group  or  of  Computer  Sciences  Department.  Other  installations  maintain  a 
library  of  computer  codes  readily  available  to  users;  NELC  slioiild  do  the  same. 

In  the  applications  area,  mathematic'il  programming  has  been  a  proved 
design  tool.  In  reference  1 ,  it  was  shown  to  be  applicable  to  typical  NELC 
problems.  Tlie  second  part  of  this  report  discussed  some  topics  whicli  were 
useful  in  obtaining  actual  solutions  to  design-relaied  MP  problems.  The  tech- 
niiiues  have  been  used  with  routines  available  at  NELC  We  have  tried  to  show 
that  MP  can  be  an  aid  to  the  design  engineer  and  not  his  replacement.  In  fact 
to  use  MP  effectively  requires  that  the  engineer  be  skillful  in  his  field  and  be 
able  to  generate  an  accurate  mathematical  modt  l  of  liis  design  problem.  We 
furtiier  recommend  that  a  short  course  or  continuing  seminar  be  offered  to 
NELC  personnel  to  familiarize  them  with  MP  techniques.  Such  a  course  was 
given  in-house  in  the  fall  of  1968  and  was  well  a'ceived. 


APPENDIX  1 :  APPLICATIONS  OF  INTEGER  PROGRAMMING 
TO  ENGINEERING  DESIGN 

In  this  appendix  we  report  on  the  resuits  of  a  literature  search  to 
detern  ine  feasibility  of  using  integer  programming  (IP)  as  a  practical  design 
aid  in  an  ongoing  NELC  task  -  BAMS  (Benchmarks  foi  Applications  of  Micro¬ 
electronics  to  Systems).  The  results  were  disappointing.  Meaningful  applica¬ 
tions  .  re  still  in  the  experimental  stage,  with  results  limited  to  relatively  small 
test  problems.  The  main  hindrance  to  successful  applications  is  not  formula- 
tional  difficulties,  but  the  lack  of  reliable  computer  codes  for  solving  the 
resulting  IP  problems.  General-purpose  IP  codes  are  severely  limited  in  the 
size  problem  they  can  solve  (a  maximum  of  60  constraints  and  60  variab.es  at 
NELC),  and  they  are  sonietimes  unreliable.  The  majority  of  the  applications 
of  IP  have  been  in  the  business  world,  and  many  algorithms  for  solving  special 
IP  problems  have  been  developed;  e.g.,  aircraft  crew  scheduling  and  warehouse 
placement.  These  methods  rely  on  the  special  structure  of  the  IP  problem 
under  consideration  and  have  worked  well.  If  the  design  engineer  is  lucky, 
his  IP  problem  may  fit  one  of  these  special  methods  (it  is  still  an  art  to  match 
the  computer  code  to  the  posed  IP  problem);  otherwise  he  must  rely  on  the 
general  IP  codes.  Here,  application  is  ahead  of  theory. 

We  give  a  briel  example  which  illustrates  the  computational  difficulties 
of  IP,  a  review  of  the  state  of  the  art  of  IP  as  related  to  BAMS,  and  finally  a 
detailed  example  in  which  IF  is  used  to  solve  a  backboard  winng  problem. 

Most  of  the  IP  work  has  been  done  in  the  linear  case;  i.e., 

Minimize  fix)  =  c  j  x  j  +  02^2  S’^n 

subjectto  g|(x)  =  a|jX|  +  a  j  2Xj2  +  • . .  +  aj^Xi  <  bj 


gm'^^^^ml'^l  +  +  •  •  •  +  ^mn^n  ^  I’m 

xj  >  0  and  Xj  an  integer 

Thus,  the  integer  linear  orogramming  (ILP)  problem  is  just  the  LP  problem 
with  the  additional  constraint  that  the  solution  be  integral.  Gomory*^^  has 
developed  an  algoritliii.  ..liich  theoretically  solves  the  above  problem  in  a 
finite  number  of  steps  N  however,  N  can  be  a  large  number  and  hence 
impractical  for  some  problems.  Since  this  area  of  mathematics  is  so  useful, 
many  methods  (some  heuristic  in  nature)  for  solving  special  ILP  problems 
have  been  reported  which  work  well.  One  metho('  which  appears  obvious, 
is  to  solve  the  associated  LP  problem,  then  round  to  the  nearest  (in  some 
sense)  integer-valued  vextor  and  use  that  for  the  solution.  Many  times  this 


will  work  fine,  if  it  is  not  too  critical  to  have  the  optimum.  If  it  is  necessary 
to  have  the  optimum  and  know  that  it  is  the  optimum,  then  other  method" 
must  be  used.  The  following  example  points  out  seme  difficulties  which 
can  arise. 

Minimize  z  =  f(Xj ,  X2)  =  -X|  -  4x2 

subject  to  0<X|<4.55 
0<X'.<4.G 

.5xj  +  X2 5.2  xj,  X2  integei-s 
The  constraint  set  (denoted  by  +)  looks  as  follows  (fig.  7): 


When  the  associated  LP  problem  is  solved  (by  inspection),  we  obtain  (1.3, 
4.55).  Rounding  to  the  nearest  vector  with  integer  components  yields  (1 ,5), 
but  (1,5)  is  not  feasible.  If  we  take  the  closest  feasible  point  to  (1.3,  4.55)  in 
the  Euclidean  norm  sense,  then  we  have  (1,4)  for  a  solution,  which  yields 
z  =  f(  1 ,4)  =  -1 7.  However,  another  feasible  point  (2,4)  yields  z  =  f(2,4)  18 

and  this  point  is  the  solution  As  the  number  of  variables  and  constraints 
increases,  it  becomes  difficult  to  find  a  feasible  point  to  the  ILP  near  the 
LP  optimum. 

The  field  of  digital  systems  design  has  made  the  most  engineering  use 
of  IP.  but  again  designers  have  had  only  limited  success.  Logic  designers  have 
used  ILP  as  a  theoretical  tool  and  have  solved  small  problems  numerically. 
Muroga'*  perhaps  has  the  most  recent  application  (see  references  52  and  53 
also).  He  discusses  designing  optimal  networks  of  the  “feedforward”  type 
by  IP  and  defines  a  generalized  gate  called  a  threshold  gate.  His  IP  formula¬ 
tion  allows  a  wide  choice  of  objective  functions,  depending  on  the  application. 


as  well  as  the  inclusion  of  any  design  constraints.  If  R  is  the  number  of  gates 
in  the  design,  then  his  associated  ILP  problem  has  variables.  Fc'  R  of  any 
typical  size  in  a  logic  design  problem,  the  resulting  ILP  problem  would  be 
intractable.  Presently,  work  in  this  area  is  of  academic  interest  only.  Also, 
with  the  abundance  of  off-the-shelf  MSI  and  LSI  components,  little  design 
is  done  at  the  gate  level  (except  by  those  manufacturing  the  aforementioned 
items).  Thus,  the  extra  effort  to  formulate  the  logic  design  problem  as  an 
ILP  does  not  appear  to  be  cost-effective. 

The  next  pertinent  application  area  is  that  of  actual  component  and 
circuit  layout.  Kodres^*^  develops  the  theory  for  solving  the  circuit  layout 
problem,  which  he  defines  as  follows; 

“The  circuit  layout  problem  is  viewed  as  a  sequence  of  four 
subproblems. 

1 .  The  determination  of  standard  replaceable  modules. 

2.  The  partitioning  of  circuits  into  groups  subject  to  input-output 
restrictions. 

J.  The  selection  of  replaceable  modules. 

4.  The  circuit  placement  and  the  interconnection  problem.” 

Kodres  uses  graph  theory,  combinatorics,  and  integer  programming  to  formu¬ 
late  the  problem.  The  actual  casting  of  parts  3  and  4  as  an  IP  problem  is 
partially  in  terms  of  some  graph  theoretical  concepts  and  requires  more  back¬ 
ground  than  we  can  present  here.  This  paper  points  the  way  for  future  work 
in  this  area;  again,  the  ideas  are  far  ahead  of  practical  methods  for  implement¬ 
ing  them.  Breuer  in  reference  55  poses  part  4  as  a  single  ILP  problem,  in 
straightforward  terms,  which  we  outline  in  the  following  paragraphs.  Other 
applications  are  in  coding  theory  and  satellite  communications  network  design; 
see  references  56  -  58. 

BREUER’S  PLACEMENT  AND  INTERCONNECTION  IP  FORMULATION 

The  backboard  wiring  problem  consists  of  three  subproblems:  the 
placement  problem,  the  connection  problem,  and  the  routing  and  installation 
problem;  each  is  dependent  on  the  other  two.  We  discuss  the  first  two  and 
pcse  them  as  a  single  ILP  problem  which,  when  solved,  will  simultaneously 
'OiVe  both  problems. 

t- .  \CEMENT  PROBLEM 

Given  B  objects,  connect  each  object  to  a  subset  of  the  remaining 
(B-1 )  objects.  The  objects  are  constrained  to  ue  on  grid  points  which  repre¬ 
sent  the  backboard  of  a  computer,  or  any  digital  system  hookup.  The  object 
is  to  place  all  the  modules  so  that  the  hcokup  wim  is  of  minimal  length. 


5.S 


CONNECTION  PROBLEM 


Given  S  fixed  objects  which  are  to  be  made  electrically  common, 
connect  the  objects  so  that  the  total  interconnection  length  is  minimal. 

INTEGER  PROGRAMMING  FORMULATION 

Given  B  objects  to  be  placed  at  the  intersections  of  the  rectangular 
grid  (fig.  8). 


X  AXIS 

Figure  8.  Back  plane  grid. 

let  Xj,  i  =  0,1 , .  .  . ,  n-1  be  the  x  coordinate  of  the  ith  object 
Yj,  i  =  0,1 , .  .  .  ,  m-1  be  the  y  coordinate  of  the  ith  object 
and  B<mn. 

The  interconnection  distance  between  the  ith  and  jth  object  is  defined 
to  be 

dij  =  |vi-Xj|  +  k|yi-yj|. 

We  note  that  no  two  objects  can  occupy  the  same  spot  at  the  same  time. 
Also  given  is  a  list  of  the  desired  connections.  We  wish  to  uniquely  position 
the  set  of  B  objects  at  the  intersections  and  determine  which  objects  should 
be  directly  connected  together,  in  a  manner  such  that  total  interconnection 
distance  of  those  objects,  directly  connected,  is  minimal.  The  hard  parts  are 
making  sure  no  two  objects  occupy  the  same  spot  and  getting  a  linear 
relation  for  d  j. 


The  constraint  that  no  two  objects  lie  at  the  same  point  requires 
that  if  Xj  =  Xj,  then  yj  +  yj  for  i  #j;or  if  yj  =  yj,  then  xj  ^  Xj  for  i 

We  list  the  required  cotistraints  for  the  placement  problem  and  then 
explain  how  they  meet  the  conditions  of  the  problem. 


Xj  -  Xj  <  n  6jj 

(Al) 

Xj  -Xj+  1  <n(l  -6jj) 

(A2) 

6jj  <  1  i.e.  6jj  =  0  or  1 

(A3) 

ajj<n5jj 

(A4) 

Xj  -  Xj<aji<Xj-Xj  +  n(l  ~5jj) 

(A5) 

aji<n(l-5jj) 

(A6) 

Xj  -  Xj  <  CKjj  <  Xj  -  Xj  +  n  5jj 

(A7) 

for  all  i  >  j,  j  =  1 ,2, .  .  .  ,  B-1 

From  inequalities  (Al)  and  (A2)  we  find  that  if  Xj  >  Xj,  then  Sy  =  1 
and  5jj  =  0,  if  Xj  >  Xj.  Inequalities  (A4)  -  (A7)  give  a  representation  for 
|Xj  -  Xj  I  as  follows.  Since  (A4)  and  (A5)  yield  0  <ajj  =  Xj  -  xj  <  n  -  1 , 
if  Xj  >  Xj  and  cvjj  =  0  otherv/ise.  In  a  similar  fashion  0  <0(j  j  =  Xj  -  Xj, 
if  Xj  >  xj;  thus,  we  have  ajj  +  Oj j  =  1  Xj  -  Xj  1 . 

Now  let  n  and  P  play  analogous  roles  for  y;  then  =  lyj  -  yj  I  • 

Thus,  djj  can  be  stated  in  terms  of  these  auxiliary  variables  as 

dij  =  aij+aji  +  lc«3ij  +  <3ji) 

Some  additional  constraints  guarantee  realizability. 


(A8) 

Xi<n-!  i-1,2 . B 

(A9) 

yj  <  ni-l  i  =  1 ,2, .  .  . ,  B 

(AlO) 

Inequality  (A8)  guarantees  that  tA'o  objects  do  not  occupy  the  same  grid  inter¬ 
section.  All  Xj  and  yj  are  nonnegative  integers.  Note  that  the  above  ipequalities 
do  not  depend  on  how  the  objects  are  connected,  but  only  guarantee  that 
specific  necessary  conditions  will  be  satisfied  Inequalities  (AI )  through  (AlO) 
represent  the  constraints  for  the  placement  problem.  We  now  turn  to  the 
problem  of  optimally  interconnecting  the  B  objects. 

Assume  that  there  are  C  independent  circuits  where  the  jth  circuit  can 
be  connected  in  Pj  different  acceptable  ways.  Let  fjj(d)  be  an  expression  for 
the  total  length  of  wire  in  the  ith  way  of  connecting  the  jth  circuit.  For 
example,  if  the  sixth  circuit  consists  of  three  objects  (1.2,  .2)  and  the  first 


may  have  only  one  connection  to  it,  then 


*^1 6"  ^21  ■^‘’•32  f26^‘*31  ■'■‘^32 

The  tv/o  final  constraints  which  relate  to  the  connection  problem  are: 

Vj^fij(d)  +  (7ij-l)Q  (All) 

and 

Pi 

^  7ij  =  l  forj  =  l...C  (A12) 

i=l 

Equation  (A12)  implies  that  each  circuit  is  connected  and  7jj  is  0  or  1. 

In  (A1 1)0=  max.  fjj(d)  fr  i  =  1, 2, . . . ,  Pj;  j  =  1, 2, . .  . ,  C.  Then  the 
resulting  objective  function  is 

C 

min.  Z  =  Vj 
J 

j=l 

The  objective  function  and  inequalities  (Al)  -  (A1 1)  form  the  ILP  repre¬ 
senting  the  combined  placement  and  connection  problem.  Although  the 
formulation  seenis  straightforward,  the  resulting  ILP  problem  can  be  large, 
even  for  small  values  of  B.  Breuer  gives  the  following  relations  between  the 
number  of  variables  I,  the  number  of  constraints  W,  and  B  for  the  placement 
problem;  i.e.,  constraints  (AD  -  (AlO). 

1(B)  =(B/2)(19B-15) 

W(B)  =  B(3B-1) 

Thus,  for  B=5  (which  could  be  manually  positioned  quickly  and  most  likely 
optimally),  we  have  1(5)  =  200,  V/(5)  =  70;  for  B  =  10  (stiil  not  too  large), 

I(  10)  =  875,  W(  10)  =  290.  These  ILP  problems  are  very  large  and  are  beyond 
the  capabilities  of  today’s  methods. 

It  seems  that,  unless  newer  formulational  techniques  are  developed 
which  lead  to  smaller  IP  problems,  or  the  computational  capabilities  solving 
IP  increase  rapidly,  this  appioach  wih  remain  a  theoretical  tool  with  no 
practical  applications.  The  following  quote  (Glover*’^,  p.  1 )  sums  up  the 
current  state  of  affairs  in  integer  programming:  “Since  its  inception  integer 
linear  programming  has.  paradoxically,  been  a  source  of  both  promise  and 
disappointment.  Promise  because  there  are  manifold  and  compelling 
opportunities  for  its  application;  disappointment  because  it  has  made  only 
the  most  dubious  proga'ss  in  spite  of  these  opportunities.” 


APPENDIX  2:  USER  INFORMATION 


HI  -  H2  OPT4LG 


CATALOG  IDENTIFICATION ; 

Ki  -  H2  0PTALG 

PROGRAMMER: 

F.  S.  Hillier,  Stanford  University,  adapted  for  NELC  by  D.  Klamer, 
Decision  and  Control  Technology  Division. 


PURPOSE: 


An  algorithm  for  solving  the  cure  integer  linear  programming  problem, 
n 


Maximize 


j-1 


subject  to 


n 


(i) 

(i  =  l,2, .. 

.  ,  m) 

j=l 

(ii) 

v 

o 

U=l,2,.. 

-  .  n) 

(iii) 

Xj  is  an  integer 

(j=i.:,.. 

.  ,  n) 

RESTRICTIONS  AND  LIMITATIONS : 

The  dimensions  of  the  constraint  matrix  A(LJ)  =  ay  have  to  bi-  less  than 
or  equal  to  61  X  61 :  i.e.,  m  <  61  and  n  <  61 . 


LANGUAGE: 

FORTRAN  IV 

COMPUTER  CONFIGURATION : 
Go  step  REGION  336k 
IBM  360/65 


METHOD: 

An  initial  noninteger  solution  must  be  obtained  from  the  related  li  near 
programming  problem  (i.e.,  xj  is  not  necessarily  an  integer)  as  well  as  tl 
resulting  basis  inverse.  To  accomplish  this,  we  have  chosen  the  linear  program¬ 
ming  routine  MPS/360.*^  Data  are  taken  directly  from  MPS/360  and  ii  put 
directly  into  Hillier's  program,'^'  without  user  intervention,  in  one  mult  step 
computer  run. 

The  advantage  of  using  this  modified  version  of  Hillier’s  program  is 
the  time  saved  from  obtaining  the  basis  inverse  and  optimal  >'.)lution;  the  time 


spent  for  punching  these  input  cards  is  also  savod.  The  data  needed  are  exactly 
the  same  as  the  first  three  card  groups  of  Hillier’s  program.  These  are: 

Card  group  1  Any  alphanumeric  characters  to  identify  the  prob¬ 
lem  in  20A4  Format 

Card  group  2  ni,  n,  KL  in  315  Format  where  A(1,J)  is  of  size  mxm 
and  KL  =  1 

'  ird  group  3  the  a  ••ays  A,  b,  c  in  15  F5.0  Format.  A(I,J)  is  the 
constraint  matrix,  B(l)  is  the  right-hand  side,  and 
C(J)  is  the  objective  function.  The  A  is  read  in  one 
row  at  a  time. 


HILLIER'S  PROGRAM 


Setup  data  for  M PS/ 3 60. 

Constraint  matrix,  right-haiui  side,  and  objective 
function. 

Using  RFADCOM  from  MPS/360  obtains  basis  inverse, 
JPM(l),  optimal  solution,  and  starting  integer  solution. 
Stores  information  on  disk. 

Computes  basis  inverse  and  opiimul  noninteger  solution. 
Computes  the  optimal  integer  solution. 

fiO 


FORTRAN  DECK  I 
DATA 

FORTRAN  DECK  2 

MPS/360 

HILLIER’S  PROGRAM 


FORTRAN  DECK  1 


This  deck  sets  up  the  daia  for  MPS/360,  since  the  format  for  MPS/360 
is  long  and  cumbersome.  The  constraint  matrix  A{I.J)  is  normalized,  as  is  the 
right-hand  side  B(l);  this  is  done  one  row  at  a  time.  The  subroutine  XPUNCH 
places  the  data  into  a  disk  file  in  proper  format  for  MPS/360,  from  which 
MPS/360  reads  the  data.  A  printout  is  given  of  the  data  that  are  placed  on 
disk.  (Note:  This  data  set  is  placed  into  a  disk  file  called  FTOIFOOI.  The 
data  are  in  normalized  form,  and,  since  MPS/360  is  designed  to  find  the  mini¬ 
mum,  the  signs  of  the  cost  coefficients  (objective  function)  are  changed.) 

The  constraint  matrix  A(I,J),  the  right-hand  side  B(l),  and  the  cost  coefficients 
C(J)  are  also  stored  on  the  disk  file  called  FT02F001 . 

FORTRAN  DECK  2  (DATAHII.L) 

This  deck  is  a  temporary  update  that  is  concatenated  onto  MPS/360, 
under  the  name  DATAKILL.  It  uses  READCOMM,^*^  which  is  a  subroutine 
designed  to  augment  MPS/360  with  procedures  written  in  FORTRAN  language. 
DATAHILL  retrieves  from  MPS/360  the  basis  inverse,  the  order  of  the  basic 
variables,  the  optimal  noninteger  solution,  and  a  starting  feasible  solution  to 
the  integer  linear  programming  problem.  These  data  are  then  added  to  the 
data  from  the  first  FORTRAN  deck  on  the  disk  file  FT02F001. 

The  starting  feasible  solution  is  a  lower  bound  on  the  value  of  the 
objective  function.  To  obtain  this  feasible  solution,  we  have  chosen  to  do  the 
following:  If  the  cost  coefficient  is  positive,  round  the  corresponding  variable 
of  the  optimal  noninteger  solution  down  to  the  next  largest  integer.  If  the 
cost  coefficient  is  negative,  round  up  to  the  next  smallest  integer.  (Note: 

The  cost  coefficients  are  p'aced  into  the  disk  file  called  FT03F001  in  the 
first  FORTRAN  deck  and  are  read  by  DATAHILL:  the  constraint  matrix  and 
right-hand  side  are  also  passed.)  This  rounded  solution  is  checked  for  feasibil¬ 
ity.  If  the  solution  is  feasible,  then  proceed  to  the  next  step.  If  the  solution 
is  not  feasible,  then  iry  a  rounding  procedure  to  satisfy  the  constraint  violated. 
(User  may  set  the  maximum  number  of  iterations  or  changes  in  the  solution.) 

In  order  for  Hillier’s  program  to  be  executed,  a  feasible  integer  solution  must 
be  found;  if  such  a  solution  is  not  found,  then  Hillier’s  program  is  skipped  and 
ail  the  data  are  punched  out  on  cards.  Tnis  information  includes  the  name, 
constraint  matrix  A(I,I),  right-hand  side  B(l).  cost  coefficients  C’(J ),  the  basis 
inverse,  the  order  of  the  ba  Ic  variables,  and  the  optimal  noninteger  solution. 
The  user  may  then  supply  his  own  starting  integer  solution  and  run  the  prob¬ 
lem  directly  from  LLOAl),  using  the  punched  data  produced  in  the  last  step. 
(See  Part  II.) 

If  a  feasible  integer  solution  is  obtained  and  the  user  wishes  to  have 
the  above  data  also,  then  ■■ither  of  iwu  metliods  may  be  use  .1.  First,  tiicrc 


(i 


arc  comment  cards  m  the  program  to  punch  out  each  of  the  groups  of  data. 

All  that  is  required  is  to  remove  the  “C”  from  the  cards  in  the  program  cor¬ 
responding  to  which  groups  are  to  be  punched  out.  The  second  method 
requires  two  changes  in  the  JCL  cards.  See  the  JCL  listing  at  the  beginning 
of  the  program. 

MPS/360 

MPS/360  is  an  IBM  supplied  application  program,  “Mathematical 
Programming  System/360.”  MPS/360  obtains  an  optimal  solution  (non¬ 
integer)  from  the  related  linear  programming  problem  and  finds  the  inverse 
of  the  basis. 

HILLIFR’S  PROGRAM 

Hillier’s  program  resides  o\i  LLOAD  (a  partitioned  data  set  on  NELC’s 
360/65  d'sk  storage)  under  the  name  OPTALG. 
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This  section  covers  the  necessary  input  for  Hillier’s  p.ogram  when  run 
without  the  adapted  program  to  generate  the  data.  The  following  is  the  neces¬ 
sary  input: 

Card  input  1  Any  alphanumeric  characters  to  identify  the  problem 
in  20A4  Format;  e.g.,  “Thompson  Number  8” 

Card  group  2  m,  n,  KL  in  315  Format  (m  <  61 ,  n  <  61 )  where  m 
is  the  number  of  rows  of  the  constraint  matrix  A(1,J1 
and  n  is  the  number  of  columns 

1  if  the  basis  inverse  is  in 
KL  =  normalized  form 
0  otherwise 

Card  group  3  This  group  of  cards  contains  the  arrays  A,  b,  c;  the 
Format  is  15F5.0.  A(l,Jl  is  the  constraint  matrix, 

B(l)  is  the  right-hand  side,  and  C(J)  is  the  row  matrix 

of  the  cost  coefficients  (or  objective  function).  The 

A  array  is  read  in  one  row  at  a  time.  (For  example. 

if  m=2,  m=16,  the  cards  would  be: 

aj  j(j=l,2  .  ...  15)  in  the  first  75  columns; 

aj  1^  m  the  first  5  columns; 

a-!  i(j=l,2, ....  15)  in  the  first  75  columns: 

in  the  first  5  columns; 
b] ,  b2  in  the  first  10  columns; 

Cj(j=i,2 . 15)  in  the  first  75  columns; 

C|^  in  the  first  5  columns.) 

Card  Group  4  This  group  of  cards  contains  the  basis  inverse  in 

6F13.5  Format.  1  he  rows  are  read  in  sequentially. 

In  the  example  I’sting, 


CO 

CO 

BB,  2 

BB,3 

•  bb,. 

BBi,7 

BBi,8 

BB-)  j  .  . 

BB24 

BB2.5 

BB2.6 

BB2  7 

BB3  ^ 

(’ard  Group  5  JPM(i),  i=l,2  .  .  .  .  in  1514  Format,  where  JPM(i)  is 
the  index  of  the  i**'  basic  variable  (including  slack 
variables)  from  the  simplex  code. 

Card  (iroup  6  The  optimal  solution  t(>  the  related  linear  programming 
problem,  x(j)  in  (>F13.5  Format. 
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Card  Group  7  A  starting  optimal  solution  to  the  Integer  Linear 
Programmii.*,  problem.  XF(j),  in  6F1.7.5  Format. 


See  the  sample  problem  following.  The  correspondence  of  cards  with 
the  above  card  groups  is  as  follows; 

Card  Group  Card  in  Sample  Problem 


1 

2 

3 

4 

5 

6 
7 


1 

2 

3-12 

13-23 

24 

25-26 

27-28 


(Note:  The  first  three  card  groups  are  the  same  as  the  three  card  groups  for 
the  adapted  version  of  Hillier’s  program.) 


SAMPLE  PROBLEM 


//PKELIM  JOB  1055369,6202044, F ,Dt5t5,1000» .DKLAMEB, 

//  MSGLEVEL=1,CLASS=L 

/  MESSAGE  OCO 

//SI  EXEC  F0RTGCLG,TIME=^1,REGI0N.G0=64K 
//FORT.SYSIN  DO  ♦ 

DOUBLE  PRECISION  i A( 611 ,BUFr ER ( 6 1) 

DIMENSION  A{61,61),  R161I,  C(6l),  NAMEJ20),  SUM(61» 

574  FORMAT! 20A4) 

READ! 5, 574)! NAME! 11,1=1,20) 

575  FORMAT! IH  ,5X,20A4) 

WRITE! 2, 574)! NAME! I), I =1,20) 

WRITE!  6, 5  75)!  NAME!  I)  ,I=l,?0) 

201  FORMAT! 315) 

REA0!5,201)  M,N,KL 
WRITE! 2,201)M,N,KL 
200  FORMAT!  15F5.0) 

WRITE! 6, 503) 

DO  183  1=1, M 

READ  !5, 2001)1 lA! J),J=1,N) 

WRITE! 2,  2001)!  IA<  J) ,J=1,N) 

WRITE! 3, 2001 )! lA! J) , J=1,N) 

CALL  CORE! BUFFER ,488) 

WRITE!  8,  2002)!  lA!  J),J=1,N) 

CALL  CORE! BUFFER, 488) 

REA0!8,2003)  ! a! I , J ) , J=1 ,N ) 

WRITE!6,500)  ! A! I , J) , J=1 ,N) 

183  CONTINUE 

WRITE!6,504) 

READ  !  5,  2001)! lA!  J) ,J=1,M) 

WRITE!  2,  2001)! lA!  J) ,J=1 ,M) 

WRITE!3,20Cli! lA! J),J=1,H) 

CALL  CORE! 80.  FER,488) 

WRITE! 8,2002)! lA! J) ,J=1,M) 

CALL  CORE! BUFFER, 488) 

READ!  8,  2003)  !  3!  J) , J=1,M) 

WRITE16, 500)18! J),J=1,M) 

WRITE! 6, 505) 

READ  1 5, 2001)1 lA! J ) ,J=1 ,N) 

WRITE! 2, 2001 )!  lA!  J),J  =  1,N) 

WRITE! 3, 2001 )! lA! J) ,J=1,N) 

CALL  C0RE!BUFFER,488) 

WRITE!8,2002)! lA! J),J=1,N) 

CALL  CORE! BUFFER, 488) 

REA0!8,20G3)  !C!J),J=1,N) 

WRITE16,500)!C! J),J=1,N) 

2001  FORMAT!  15A5) 

2002  FORMAT16JA5) 

2003  FORMAT! 61F5.0) 

500  FORMAT!  IH  ,  15F8.  2) 

501  FORMAT!  IH  ,  15F8.4) 

555  FORMAT!  IH  ,  lOF 12.51 

567  FORMAT!  IH  ,  15Ffl.4) 
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503  FO<^^^AT(23H  CONSTRAINT  MATRIX  A(I,J)  IS> 

504  FORMAT!  24H  RIGHT  HANO  SIDE  BID  IS) 

505  FORMAT! 27H  COST  COEFFICIENTS  C{J)  ARE) 

506  FnRMAT(43H  THE  NORMALIZED  CONSTRAINT  MATRIX  A{I,J)  IS) 

507  FORMAT!  3RH  THE  NORMALIZED  RIGHT  HAND  SIDE  BCD  IS) 

K  =  0 

DO  111!  I=).,M 
K=X  +  1 

SUM!K)=C  J 
DO  2222  J=  l,N 
S'JM!K)  =  SUM!K)*A!  I,J)**2 
2222  CONTINUE 

SUMIK  )  =  SQRT!  SUM!  K)  ) 

1111  CONTINUE 

DO  3333  1=1, M 
DO  4444  J= 1,N 
A!  I, J  )  =  A!  I , J )/SUM(  I ) 

4444  CONTINUE 
3333  CONTINUE 

WR  ITE!  6,  506) 

DO  5555  1=  1, M 

5555  WR1TE!6,567)!A! I,J),J=1,N) 

DO  7777  1=1, M 
7777  b!  I )  =  B!  I )/3UM!  I  ) 

WRITE! 6, 507) 

WRITE!6,501)!8! I),I=1,M) 

508  FORMAT! • 1  THE  FOLLOWING  DATA  HAS  BEEN  STORED  CN  DISK  •> 
WRITE  !6,5Ce) 

C4LL  XPUNCH  !A,B,C,M,N) 

STOP 

END 


SUBROUTINE  XPUNCH  !A,B,C,M,N) 
DIMENSION  A!6L,61),  B(61),  C!6U 
WRITE  (1,111) 

PRINT  110 
WRITE  !  1, 1  13) 

PRINT  112 
WRITE  !  ),  115) 

PRiNT  115 
DO  I  KR'I.,  =  1,  M 
KRCWIO  =  KROW  +  10 
WR  ITF  (1,107)  KROWIO 
PRINT  1C6, KROWIO 
1  CONTINUE 

WRITE  I  1,  1C<:) 

PRCNT  ICfl 


Ml  =  M+ 1 

DO  10  J  =  1,N 

K  =11 

I  =  C 
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5  CONTINUE 


K1 

=  K 

K2 

=  K+1 

L 

=  J  +  10 

I 

=  I  +  1 

IF  (K2 

!  .GE.  Ml  +  10)  GO  TO  11 

WRITE 

(1,100)  L,K1,A( I , J) ,K2 ,A(I+l ,J) 

PRINT 

100,L,K1,A( I ,J)  ,K2,A(I  +  1,J) 

K 

=  K  +  2 

I 

=  I  +  1 

GO  TO 

5 

11 

IF  (K 

.EQ.  Ml  +10)  GO  TO  12 

C(  J  ) 

=  -  C( J) 

WRITE 

(  1,101)  L,K1,A( I ,J) ,C( J) 

PRINT 

101,L,K1,A( I , J) ,C( J) 

GO  TO 

10 

12 

C(J  ) 

=  -  C( Jl 

WRITE 

(  1, 102)  L,C( J) 

PR  INT 

in2,L,C( J) 

10 

CONTINUE 

WRITE 

(  1,  105) 

PR  INT 

lOA 

DO  20 

I  =  1,M 

IR  =1  +  10 

WRITE  ( 1, 103)  IR  ,B(  I  ) 
PRINT  103,  Ik  ,8(  I  ) 


20 

CONTINUE 

WRITE  (  1,117) 

PRINT  116 

100 

FURMAT(4X,‘C',I2,7X, 

•R* ,I2,7X,F12o5,3X,'R' ,I2,7X,F12.5) 

101 

PORMAT( AX, 'C • , I2,7X, 

•R',I2,7X,F12.5,3X,*C',9X,F12.5) 

102 

FORMAT(AX,»C',I2tyX, 

•O' ,9X,F12.5) 

103 

FORMAT!  AX,  •CONST'  ,'jX 

,'R'  ,I2,7X,F12.5) 

lOA 

FORMAT  (  •  RHS'  ) 

105 

FORMAT  (  'RHS'  ) 

106 

FORMAT!  •  L  R‘,I2) 

107 

FORMAT  (  •  1.  R',I2* 

108 

FORMAT  (  '  COLUMNS') 

109 

FORMAT  (  'COLUMNS'  1 

110 

FORMAT  ('  NAME 

0PT8AS') 

111 

FORMAT  (  'NAME 

OPTBAS') 

112 

FORMAT  (•  ROWS') 

113 

FORMAT  (  'ROWS') 

HA 

FORMAT! •  N  O') 

115 

FORMAT!  •  N  O'J 

116 

FORMAT  (•  ENOATA') 

117 

FORMAT  (  'ENOATA') 
RETURN 

END 

67 


//GO.FTOlFOOl  Dl)  UN  I  T=SYSnA  ,  SPACE=(  TRK,  (20 ,101  ,RLSE)  , 
//  DISP=( ,PASS), 

//  0Ctl^(RECFM  =  FBS,Lr<ECL=3C,8L:<SI  ZE=880i 

//GO .FT02FCni  on  UNIT=SYSDA, SPACE=(TRK, (20,10) ) , 

//  DISP=(  ,PASS)  , 

//  0C3=  (RECFM  =  F0,LRECL=^8C,BLKSIZE  =  880) 

//GO  .FTC3F001  OD  UN I  T=  SY  SO  A  ,  SPACE  =  (  TRK  ,  (  20 , 1 C  )  ,RLf.  E  )  v 
//  l)ISP=(  ,PASS)  , 

//  nCB=  (kECE'1=FB'3,LRECL  =  30,BLKSI  2E=880) 

//GO.  SYS  IN  DO  * 


$!(.$$$ 

i 

i  DATA 
S 


//STEPQNE  EXEC  FURTGCL 
//FORT.SYSL  IN  DQ  D  I  SP--(  NE  W ,  PA  SS) 

//FORT,  SYS  IN  DD  * 

INTEGER  OIFFfFILE 
INTEGER*2  CQMP( 61) ,C,R 

DIMENSION  BASE! 61,61) ,NROW( 61)  ,NC0LM(61)  ,  JPM(6l ) , I  ROW  (61 ) , 
I  BUFFER! bl),CC( 61) ,XF(6l) 

DIMENSION  A! 61,61) ,a( 61)  , 

1  IDXU  61 )  ,  IDXF!  61)  ,IDX2(  61  )  ,  XFMIN(61  ) 

REAL^a  NAME,XLIST(  30  ,  OUTl  (  62  )  ,  0UT2  ( 372  ]  ‘  ,  OUT  16  1 ) 

DATA  BASE/372i*0,0/ 

LK  =  0 
F  ILE=4 
C 

C  TO  OBTAIN  THE  NUMBER  OF  ROWS  AND  COLUMNS  CF  THE 
C  CONSTRAINT  MA  TR  I  X(  I  NCLUDI  NG  THE  OBJECTIVE  FUNCTION) 

C 

CALL  ARRAYiFlLE,INniC,NAME) 

CALL  vECTORIFILE, INDIC ,XLI S7) 

MROWS  =IFIX(  SNGL!  XLI  ST»9)  M 
NCOLMN= IF IX( SNGL!  XLI ST( 10) ) ) 

MROWSl  =  MROwS-  1 
MN1=MR0WS1  *  NCGLMN 
M;M  =  NCOLMN*MROwS 
C 

C  TO  OBTAIN  THE  BASIS 
C 

CALL  ARRAYIFILE,  INniCfNAME) 

DO  10  J=1,MN 

call  VECTOR! F ILE , INDIC , XLI ST) 

IF! J ,GT. MROWS)  GOTO  12 
O  JTl!  J  )  -  XL  I  ST!  1  ) 

12  1F!M0D! J. MROWS  )-l)  11  ,115,11 
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11 


LK  =LK  I 
0UT2(LK)  =  XLIST(2) 

115  IF(INOIC-l)  2C,2C,10 
10  CONTINUE 

20  CONI  INUE 

C 

C  TO  OBTAIN  THE  OPTIMAl  SOLUTION 

CALL  ARRAY! F ILE , INOIC , NAME) 

CALL  VECTOR(FILE,INOIC,XLI  ST) 

CALL  ARRAY(FILE,INDICfNAME) 

CALL  VECTOR(FILEiINOIC,XLIST) 

CALL  ARRAY! FILE, INOIC , NAME) 

DO  201  J=1,NC0LMN 

CALL  VECTOR!FILE, INOIC, XLI  SI) 

OUT!  J  )^^XL!  ST!  3) 

IF!INDIC-1)  200,200,201 
201  CONTINUE 
C 

C  TO  PRINT  OUT  THE  OPTIMAL  SOLUTION 
C 

101  FORMAT! 6F13. 5) 

102  FORMAT!  IX, 6F13,5) 

103  FORMAT!'  THE  OPTIMAL  NONINTEGER  SOLUTION  IS  •//) 

WRITE! 6, 103) 

WRITE!  6,  10  2)  !0'JT!  J)  ,J  =  l  ,NCOLMN) 

C 

C  TO  DETERMINE  WHICH  VARIABLES  ARF  ACTIVE 
r 

CALL  CORE! BUFFER, 240) 

WRITE!8,104)  ! OUT  1 !J ) , J=2 , MRGWS) 

104  FORMAT161A4) 

CALL  CORE!  BUFFER, 2‘-C) 

READ  !8  ,105)  ! COMP ! J ) , I  ROW ! J ) , J  =  1 . MROWSl ) 

105  F0RMAT!61!  1A1,I2,  IX)  ) 

DATA  R/'R  •/,C/'C  '/ 

NiiROW  =  0 

NQCOLM  =  0 
DO  500  J=1,MRQWS1 

IF!COMP! J)  .FQ.R)  NOROW  =  NOROW  ♦  1 
500  IF!COMP! J  )  .EQ.C )  NOCOLM  =  NOCOLM  +  1 

IF! !NORUW.EU.  0  ) . A  NO . ! NOC OLM. E U.  0  ))  GOTO  100 
IF!NOROW.Eg.  0  )  GOTO  600 
IF!NOCQLM.Eg.  0)  GOTO  650 
C 

C  TO  COMPUTE  THE  J TH  BASIC  VARIABLE  FOR  OPTALG  USING  MPS/360 
INFORMAT  I 
C 

DO  56  J=l, NOROW 
56  NROWIJ)  =  IROWIJ) 

I  =  0 

NTOT  =  NOROW  NOCOLM 
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Ni  =  NUROW  ♦  1 
DO  57  J=Nl,NTOT 
I  -  1+  I 

57  NCOLM(  I  )  =  IROW( J) 

GOTO  150 
C 

C  BASIS  INVERSE  IS  BASI  S 
C 

600  00  601  J=l,NOCOLM 

601  JPM(J)  =  J 
WR ITE( 6, 105) 

DO  602  I=l,MRawSl 

602  WRITE(6,101)  (  OUT2I  I  ♦  t  J-1 1  ♦HROWSl)  t  .?=1  t  NCOLPNI 

WRITE! 2,101)  ( (OUT2( !♦ { J-1 ) ♦MROWSl ) , J=1 , NCOLPN) , 1= 1 , MRCWS 1 ) 
C 

C  IF  YOU  WANT  THE  BASIS  INVERSE  PLNCHED  OUT 
C  (THERE  ARE  THREE  (3)  CARDS  FOR  THE  BASIS  INVERSE) 

C  REMOVE  THE  'C  FROM  THE  FOLLOWING  CARD* 

C  WRITE(7,101)  ({0UT2{I«-(J-1)*MR0WS1),J=1,NC0LPN),I  =  1,MRCWSI) 
C 

GOTO  202 
C 

C  BASIS  INVERSE  IS  IDENTITY  MATRIX 
C 

650  DO  651  J=l, NOROW 
JPM(J)  =  MROWSl  *■  J 

651  BAS£(J,J)  =1.0 
WRITE! 6, 105) 

WRITE (6, 101)  (!BASE!I,J) ,J=1, MROWSl) ,1=1, MROWSl) 

WRITE (2, 101)  (! BASE! 1, J) ,J=1, MROWSl)  ,1=1, MROWSl) 

C 

C  IF  YOU  WANT  THE  BASIS  INVERSE  PLNCHED  OUT 
C  (THERE  ARE  THREE  (3)  CARDS  FOR  THE  BASIS  INVERSE) 

C  REMOVE  THE  'C*  FROM  THE  FOLLOWING  CARD. 

C  WRITE(7,101)  ( (BASE! I ,J) ,J=1 , MROWSl)  ,1=1  , MROWSl) 

C 

GOTO  202 

150  CONTINUE 
L=1 
KL=  1 

NORMl  =  NOROW  -  1 
IFINROW!  D-  11)  100,30,22 
30  JPMIl)  =  NROWII)  ♦  NCOLMN  -  10 
L  =  Ltl 
GOTO  57 

22  DIFF  =  NROW!  1)  -  11 

DO  46  K=1,0IFF 
JPMIL  )  =  NCOLM(KL)  -  10 
KL  =  Kl*l 
L  =  L*1 

46  CONTINUE 

JPM(L)  =  NROW(l)  ♦  NCOLMN  -  10 

L  =  L  ♦  I 

47  CONTINUE 

DO  50  1=1, NORMl 
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DIFF  =  NROW(H-l)  -  (NRQWnj+l) 

IF  (DIFF)  100,35,40 

35  JPM(L)  =  NR0W(I+1)  ♦  NCOLMN  -  10 
L  =  L  +  1 
GOTO  50 

40  DO  45  K  =  1,DIFF 

JPM(L)  =  NCOLMIKL)  -  10 
KL  =  KL+1 
L  =  L  +  1 
45  CONTINUE 

JPM(L)  =  NR0W(I  +  1)  +  NCOLMN  -  10 
L  =  L  +  l 
50  CONTINUE 

DIFF  =  10  ♦  NTOT  -  NROWtNOROk,) 

IF(DIFF)  100,55,60 

5^  JPM(NTOT)  NR0W(  NOROW)  *■  NCOLMN  -  10 

GOTO  65 

60  DO  65  K=l,DIFF 

JPM(L)  =  NCOLM(KL)  -  10 
KL  =  KL+1 
L  =  L  +  l 

65  CONTINUE 

C 

C  TO  OBTAIN  THE  BASIS  INVERSE 
C 

J  1  =  0 
J2  =  1 

DO  300  JK=l,NTnr 

IF( JPM( JK ) .GT. NCOLMN)  BASE(J2,JK)  =  1,0 
IF( JPM(JK) .LE. NCOLMN)  GOTO  310 
J2  =  J2  +  1 
GOTO  300 

310  DO  311  J^UMROWSl 
J5  =  Jl  'r  1 

311  8A:-E(  J,JK  )  =  0UT2(  Jl) 

300  CONTINUE 

C 

C  THE  FOLLOWING  IS  FOR  OUTPUT. 

C 

WRITE(  6v  109) 

WRITE!  6, 101)  ( (8ASE( I , J) , J=l ,MR0WS1)  ,I=1,MR0WS1) 
WRITE(  2,  101)  ( ( 8ASE( I , J)  ,  J  =  l ,MRnWSl)  .1=1  ,MR0WS1) 

C 

C  IF  YOU  WANT  THE  BASIS  INVERSE  PUNCHED  OUT 
C  (THERE  ARE  THREE  (3J  CAROS  FOR  THE  BASIS  INVERSE) 

C  REMOVE  THE  'C •  FROM  THE  FOLLCWING  CARD, 

C  WRITE(7,101)  (  (BASE!  I , J)  ,J  =  l ,MR0WS1)  ,1=1 .MROWSl) 

C 

202  CONTINUE 

WRITE! 6,  107) 

WRITE(6,108)  (JPMI.fK),  JK  =  l,NTOT) 

WRITE(2,108)  (JPM(JK),  JK=1,NT0T) 

C 


C  IF  YOU  WANT  THE  ORDER  OF  THE  I TH  BASIC  VARIABLEf  JPf'd)* 
PUNCHED  OUT 

L  REMOVE  THE  'C  *  FROM  THE  FOLLOWING  CARO. 

C  WRITE(7,108)  (JPM(JK),  JK=l,NTOT) 

C 

WRITE(2fl01)  (OUT( J) ,J=l,NCOLMN) 

C 

C  IF  YOU  WANT  THE  OPTIMA-  '  >j\TEGER  SOLUTION  PUNCHED  CUT 
C  REMOVE  THE  'O'  FROM  THE  FOLLOWING  CARD. 

C  WRITE(7,10n  (OUT!  J>  ,J=1  ,NCOLMN) 

C 

C 

C  IF  YOU  WANT  TO  SUPPLY  YOUR  OWN  FEASIBLE  INTEGER  SOLUTION 
C  REMOVE  THE  'C  FROM  THE  FOLLOWING  CARD. 

C  GOTO  151 

C 

C 

C  TO  OBTAIN  THE  STARTING  INTEGER  SOLUTION  TG  THE  ILP 
C 

DO  no  I=l,MROWSl 
READ!  3,111)  I  Al  I  ,  Jl  ,  J=1  ,NCOLMN) 
no  CONTINUE 

READ!  3,111)  (B( I ) ,I=l,MROWSl) 

REAO(3, 111)(CC( J),J=l,NCOLMN) 

Ill  FORMAT!  15F5.C) 

DO  1000  I=1,NC0LMN 
XF(  I)  =  SNGLIOUTdl) 

IFICCIJ))  1001,1000,1000 
1001  IF  (  XF(  1)  .EQoAINT!  XFCin  )  GOTO  1000 
XFd  )  =  XF(  I  )  +  1. 

1000  CONTINUE 

DO  1020  I=1,NCULMN 

1020  XF(  I  )  =  AINT!  XFd  )  ) 

ITEST  =  0 

C 

C  SET  THE  MAXIMUM  NUMBER  OF  ITERATIONS. 

C 

1021  IFCIO  -ITEST)  1100,1022,1022 

1022  ITEST  =  ITEST  +  1 
C 

C  CHECK  TO  SEE  IF  THE  POINT  IS  FEASIBLE. 

C 

DO  1030  I=l,MROWSi 
SUM  =  0.0 

DO  1025  J^1,NCULMN 

SUM  =  SUM  ♦  A! I , J)  ♦  XF( J) 

1025  CONTINUE 

IF(Bd)-SUM)  1026,1030,1030 

1026  ITER  =  I 
C 

C  IF  POINT  IS  NOT  FEASIBLE,  WRITE  THE  CONSTRAINT  VIOLATED  AND 
C  THE  POINT  THAT  VIOLATED  THE  CONSTRAINT. 


c 

WRITE(  &♦  1027)  I,B(  1)  ,SUM,(  XF(IJ)  ,IJ=1  jNCCD'N) 

1027  FORMAT!  10X,10H  B(I2t21H)  SUM  / 

1  15X,  2ri2.2  //  14H  THE  POINT  IS  /4 (15 F8 .2/ )////// ) 

GOTO  10  35 
1030  CONTINUE 
GOTO  1120 

1035  L  =  0 
N1  =  0 

C 

C  FIND  ALL  OF  THE  NONZERO  VALUES  OF  THE  VARIABLES. 

C  IDXl  STORES  ALL  OF  THE  ZERO  VALUED  VARIABLES. 

C  IDXF  STORES  ALL  OF  THE  NQNZ'RO  VALUED  VARIABLES. 

C 

DO  1040  J=1,NC0LMN 

IF!  XF!J)  )  1036,1037,1038 

1036  XF!J)  =  0.0 

1037  L  =  L+1 
IDXKL  )  =  J 
GOTO  1040 

1038  N1  =  N1  ♦  1 
XFMIN!N1)  =  XF!J) 

.  IDXriNl)  =  J 

1040  CONTINUE 
IF!L.LT.l)  GOTO  1042 

C 

C  TACK  ON  TO  THE  END  OF  THE  INDES  OF  THE  NONZERO  VARIABLES  THE 
C  INDEX  OF  THE  ZERO  VARIABLES. 

C 

DO  1041  1=1, L 

1041  IDXF!N1  +  I  )  =  IDXl!  1) 

C 

C  IF  THERE  IS  AT  LEAST  ONE  NONZERC  VARIBALE, 

C  THEN  ARRANGE  THE  NONZERO  VARIABLES  FROM 
C  MINIMUM  TO  MAXIMUM  VALUE. 

C  inX2  IS  the  INDEX  OF  THE  REARRANGED  VARIABLES. 

C 

1042  IFINl.GT.O)  GOTO  1045 

DO  1044  1=1, L 

10A4  XFMIN!n  =  0.0 
GOTO  1046 

1045  CALL  SORTI! XFMIN, 10X2,1, Nil 

1046  IF!N1.EQ.NCULMN)G0T0  1048 
NlPl  =  N1  ♦  1 

DO  1047  I  =  NIP1,NC0LMN 

1047  10X2!  I)  =  IDXF!  I  ) 

C 

C  ROUNDING  PROCEDURE. 

C  TAKE  THE  SMALLEST  NONZERO  VARIABLE, 

C  IF  THE  CORRESPONDING  CONSTRAINT  COEFFICIENTS  IS 

C  POSITIVE  —  ROUND  DOWN 

C  NEGATIVE  —  ROUND  UR. 

C 


1048 

no  1060  KL=  1,NC0‘  MN 

IV  =  IDXF!  IDX2IKL)  ) 

LCK  =  0 

1049 

IF!  A!  ITER, IV)  )  1050,1060,1051 

1050 

XF!  IV)  =  XF!  IV)  +  1.0 
GOTO  1055 

1051 

XF! IV)  =  XF! IV)  -  1.0 

IF! XF! IV) .L T.0.0)  XF!IV) 

=  0.0 

1055 

SUM  =  0.0 

DO  1056  I=1,NC0LMN 

SUM  =  SUM  A(  ITER,I  )  ♦ 

XFII  ) 

1056 

CONTINUE 

c 

c 

C  IF  THIS  DOES  NOT  STAISFY  THE  CGNSTKAINT, 

C  THEN  TAKE  THE  NEXT  SMALLEST  VARIABLE 

C  AND  REPEAT  THE  PROCEDURE. 

C 

C  IF  THE  CONSTRAINT  IS  SATISFIED,  THEN  USE  THIS  PCINT  TO 
C  CHECK  ALL  OF  THE  OTHER  CONSTRAINTS. 

C 

C 

IF(  B(ITER)  -  SUM)  10  57,1021,1021 
1057  LCK  =  LCK  ♦  1 
C 

C  SET  THE  MAXIMUM  NUMBER  OF  CHANGES  FOR  ONE  VARIABLE  FCR  ONE 
ITERATION. 

C  CHANGE  THIS  IF  THERE  IS  AN  OSCILATION  BACK  AND  FORTH  BETWEEN 
C  TWO  POINTS. 

C 

IF(LCK.LT.l)  GOrO  104? 

1060  CONTINUE 
GOTO  1021 
1100  CONTINUE 

WRITE!  6,  1105)  ITEST 

1105  FORMAT! 47H  NO  FEASIBLE  STARTING  INTEGER  SOLUTION  HAS  BEEN 
1  26H  FOUND  AT  THIS  POINT  AFTER  15, 6H  TRYS.) 

GOTO  1150 
1120  CONTINUE 
OBJ  =  0.0 
DO  1122  I=l,NCOLMN 
1122  OBJ  =  OBJ  ♦  XF!I)*CC!I) 

WRITE! 6,  1125)  ! XF! 1 ) ,1=1 ,NCOLMN) 

1125  FORMAT!53H  A  FEASIBLE  STARTING  INTEGER  SOLUTION  HAS  BEEN 

FOUND.  / 

14! 15F8.^/)  ) 

WRITE!6,  1126)  OBJ 

1126  FORMAT!  40H  THE  VALUE  OF  THE  OBJECTIVE  FUNCTION  IS 
/  F12.4  ) 

WRITE!2,101)  !  XF!  J)  ,  J  =  l,NCOLMN) 

C 

C  IF  YOU  WANT  THE  OPTIMAL  'FEASIBLE*  INTEGER  SOLUTION  PUNCHED  OUT 
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C  REMOVE  THE  'C*  FROM  THF  FOLLOWING  CARD. 

C  WRITE(7,10U  (  XF(  J)  ,J=l  ,NCOLMN) 

C 

GOTO  151 

100  WRITE(6,106J 

106  FORMAT!  *  AN  ERROR  HAS  QCCUREOM 

107  FORMAT! •  JPM! 1)  IN  THE  ORDER  THEY  SHOULD  OCCUR  AND  IS  THE 
INDEX  OF  1  the  ITH  BASIC  VARIABLE*) 

108  FORMAT! 151^) 

109  FORMAT!////*  THE  BASIS  INVERSE  IS  *//) 

1150  CONTINUE 
I  =  10 

WRITE! 10,1151)  I 

1151  FORMAT! 15) 

151  CONTINUE 

I  =  0 

WRITE!  10,1)51)  I 
200  RETURN 
END 

SUBROUTINE  SORTI ! A , 1 0 , 1  I , J J) 

DIMENSION  A! 1) ,IU!  30) ,IL! 30) ,I0!1) 

INTEGER  T1,T2 
M=1 
I=II 
J  =  JJ 

DO  1  IS=I,J 
1  ID!IS)=IS 

5  IF!  I  .GE.  J  )  GOTO  70 

10  K=I 

IJ=! J+I)/2 
V=A(  IJ  ) 

Tl=  I0(  IJ  ) 

IF!A!  I)  .LE.  T)  GOTO  20 
A!  IJI=A!  I) 

ID!  IJ)=ID(  I) 

A!  I  )=T 
ID!  I  )=T1 
T=A!  IJ  ) 

T1=ID!  IJ) 

20  L  =  J 

IF!A!J)  .GE.  T)  GOTO 
A! IJ  )  =  A! J  ) 

ID!  IJ)=ID!  J  ) 

A!J)  =T 
ID!  J  )=T1 
T  =  A(  IJ  ) 

Tl=  ID!  IJ  ) 

IF(A!  I)  .LE.  T)  GOTO  40 
A!  IJ  )  =  A!  I  ) 

ID!  IJI=IO!  I  ) 

A! »  )  =  T 
ID!  I  )=T1 
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T  =  A(  IJ  ) 

Tl=  ID(  IJ  ) 

GOTO  40 

30  A(L)  =  A(K) 

ID(L )=  ID(K  ) 

AIK  )  =  TT 
IDIK 

40  L=L-1 

IF! AIL )  .GT.  TI  GOTO  40 
TT=A(L  ) 

T2=  IDIL  ) 

50  K=K+1 

IFIAIK)  .LT.  T)  GUTO  50 
IFIK  .LE.  L)  GOTO  30 
IF(L-I  .LE.  J-KI  GOTO  60 
ILIM  )=  I 
lUIM  }  =  L 
I=K 

M=M  +  1 
GOTO  80 
60  IL(M)  =  K 

lUIM 
J  =  L 
M=M+1 
GOTO  80 
70  M=M-1 

!F(M  .EQ.  0)  KETURM 

I=IL(M) 

J=  IU(M) 

80  ir-<J-I  .GE.  11)  GOTO  10 

IF(  I  .EQ.  II)  GOTO  5 
1=  I-l 
90  1=1  +  1 

IFI  I  .EQ.  J)  GOTO  70 
T=A(  I  +  l) 

Tl=  ID!  I  +  l) 

IFI A(  I  )  .  LE.  T)  GOTO  90 
K=  I 

IDIK  )=  lOI  I  ) 

100  AIK+1)=A(K) 

IDIK  +  1)=  IDIK) 

K=K-  1 

IFIT  .LT.  AIK))  GOTO  ICO 

AIK  +  1)  =  T 

IDIK  +  l)  =  Tl 

GOTO  90 

END 


//LKEO.SYSL  IB  OD  0 SN  =  LLOAD ,0  I SP= I SHR ,KEE P) 

//  00  OSN=SYS1.FORTLIB,DI SP=SHR 

//LKEO.SYSLMOO  00  OSN  =  «.MCCAL  tOI  SP=I  NE  W  ,  PASS)  tUNI  T  =  SYSDA, 
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It  SPACE=(CYL,( I, If  10) ) ,DCB=(DS0RG=P0fRECFM=U,BLKSlZE=3625 ) 
//LKE0.SY.'>IN  00  ♦ 

INSER'.  READCOMM 
ENTRY  MAIN 
NAME  DATAHILL(R) 


//CPC  EXEC  PGM=COMPILER 

//STEPLIB  OD  DISP=( SHR.KEEP) tDSNAME=LLCAO 

//SCRATCHl  DO  UN  I  T=  SYSOA  t  f)  I  SP--(  NEW  t  OELF  TE)  *  S  PACE- ( TRK  i1  ]ii 

//SCRATCH2  DO  UN  IT=SYSOA  .01  SP  =  (  NEwloELETE)  :SP^CE=(  J5lS' !  ' 

//SCRATCH3  00  UN  IT- SYSOA  .0 1 SP  =  ( NEW .DELE TE ) , SP ACE= ( TR K* ( 1  *  1  I  J 

//SCRATCH4  DO  UN  I T=SYS0A ,0  I SP -( NEW .Of LFTE) , SPACE= (TRK* ( 1  *  1 )  ) 

//SYSPUNCH  00  SYSOUT  d 
//SYSIN  00  ♦ 

PROGRAM  ( 'ND* ) 

INITIAL/ 

MOVE! XOATA , 'OPTaAS' ) 

MOVE! XP8NAME, *MYF iLE ') 

MOVE! XOBJ, 'O'  i 
MOVE!  XRHS.  'CONST* ) 

ASSIGN  (  'COMMFMT*  ,*FT04F001*,*CCNM*) 

PREPOUTI  ‘CQMMFMT*  ) 

CONVERT! ' SUMMARY* ) 

BCnOUT 
SETUP!  I ) 

PR  IMAL 
SOLUTION 

TRANCGL  !  '  ENTIRE  *  » •  I  N  VER  SE  * ) 

TRANCOL!  'ENTIRE  * ) 

TRANCOL  !  'F  ILE  * , 'COMMFMT* , 'ENTIRE •) 

SOLUTION!  'FILE'.'COMMFMT*) 

DATAHILL 

EXIT 

PEND 


//EXEC  EXEC  PGM=EXECUTOR,CaNO=!C,NE,CPC) ,REGI0N=220K,T IME=3 
//STEPLIB  00  n ISP=! SHR.KEEP) ,OSNAME=LLCAO 
//  00  OSN^CMCCAL.OI SP^IOLOfPASS) ,UNI T=SYS0A 

//SCRATCHl  00  UNIT=SYS0A, OISP  =  !NFW, DELETE)  ,SPACE  =  (CYL, (1,1)  ) 
//SCRATCH2  00  UN  I T= S YSOA , 0  I SP  =  ( NE W ,DE LE T6 )  , S PACE^ ( CY L , ( 1 , 1 )  ) 
//PROBFILE  00  UNIT=SYS0A,SPACE  =  (CYL,(1  ,1)  )  ,L)!SP=(NEW, DELETE) 
//ETAl  00  UNIT=SYSOA, SOACE  =  (CYt., (1, 1)  )  ,OISP=(NEW, DELETE) 

//MATRIXl  OD  UNIT=SYSOA,SPACE  =  (CW,(  1  ,1)  )  ,l)ISP=(NEW, DELETE) 
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//SYS'^LCP  DO  UNlT=SYSOA, OSNAME=*. CPC.  SYSMLCPfDlSP=(OLO, DELETE) 
//SYSPHINT  on  SvSOUT=A 

//EXtC.FT02r0Cl  DD  D . S l.GO.F T0?F 001 ,01 SP= ( POO  i  F ASS ) 

//EXEC  .FT03F001  DD  0 SN=^ . SI .GO. F T03F 001 ,01 SP= < CLO t DELET E ) 
//FXEC.FTObFOOl  DO  SYSOUT=A 
//EXEC.FT07FOOI  DO  SYSOUT=B 

z/eXEC .FTOAFOOl  OD  UN  I T= SYSD A , SPAC E= (C YL  ,  ( 1  , 1 1 ) 

//SYSPUNCH  DO  SYSnUT=B 

//eXEC.FTlOFOOl  Of)  UNIT  =  SYSOA,SPACE=(TRK,(l,l)  ,RLSE)  , 

//  DISP=l ,PASS>, 

//  DCB=  (RECFM=FtiS,LRECL=80,BLKSI  ZE=880) 

//SYSIN  DO  DSN  =  *. SI. GO. FTOlFOOl, DI  SP=(OLO, DELETE) 

//CHKF  EXEC  FURTGCLG, REGION. G0=42K, TI ME=1 
//FORT  .SYSL  IN  UiJ  U  I  SP  =  I  NEW.PASS) 

//FORT, SYSIN  DO  * 

READ!  10,100  I 
100  FORMAT! 15) 

IF! I )  20, 20,10 
10  STOP  10 

20  CONTINUE 

STOP 
END 

//GO.FTIOFOOI  DU  D  SN=*  .E-XEC  .  F  T1  OF  001  ,CI  SP=  ( OLD  ,DE  LET  E) 


//TWO  EXEC  PGM=nPTALG,COND  =  12,LT,CHKF.GO)  ,REGION=(  ,336K),TIHE=2 
//STEPLIB  DO  OSN=LLOAO,DI  SP=SHR 
//FTCGFCOl  DO  SYSOUT=A 

//FT05F00)  DO  f) SN  =  * .  S 1  .GO.F  T02F  001 , D I  SP=!  OLD ,DE LETE ) 


//SFAIL  EXEC  rORTGCLG,COND.FORT=!0,EO,CHKF.GO) , 

//  CCINO.LKEO=!  1  0,  EQ,CHKF  .GO)  ,  I  A  ,  LT  ,F  ORT )  )  , 

//  CUND.GO=!  !  0,  EQ,CHKF  .GO  »  ,  (  4  ,  LT  ,F  OR  T)  ,  !4,LT,LKED)  )  , 

//  REGION. G0=42K,T1ME^I 
//EORT.SYSL  IN  CO  0  I  SP=!  NE  V,,PASS) 

//FORT. SYSIN  DO  ^ 

REAL  X!8C) 

5  READ! 5, iO,END=lOO)  X 

10  FORMAT!80An 

WP^ITE!7,l0)  X 
GOTO  5 
too  STOP 
END 

//GO.FTOSFOOl  DD  0 SN  =  * . S 1 .GO.F T02F 001 ,i  I >P* ( OLD , DE LETE ) 


NELC  ZANGWL 


CATALOG  IDENTIFICATION: 

E4  NELC  ZANGWL 

PROGRAMMERS: 

D.  C.  McCALL,  Decision  and  Control  Technology  Division,  and 
C.  M.  BECKER,  Applications  Software  Division 

PURPOSE: 

To  compute  the  minimum  of  a  function  f(X] , . . .  ,  x^^),  of  n  real  variables 

RESTRICTIONS  AND  LIMITATIONS: 

A  maximum  of  20  variables  can  be  handled. 

LANGUAGE: 

FORTRAN  IV 

COMPUTER  CONFIGURATION: 

IBM  360/65 

Core  storage:  19086  bytes 

ENTRY  POINTS 
ZANGWL 

SUBPROGRAMS  AND  WHERE  REFERENCED : 

User-supplied  programs 

FUNC  called  by  ZANGWL,  (POWELL) 

Programmer-supplied  programs 
ZGITER  called  by  ZANGWL 
POWELL  called  by  ZANGWL 

USAGE: 

CALL  ZANGWL  (XI,  N,  EACCUR,  QSTEP,  ISTOP,  LPRINT.  IX, 
LPUNCH,  XOPT,  FF) 

For  a  description  of  parameters  see  the  listing. 

INPUT  FORMAT: 

All  input  is  through  the  parameter  list  except  when  user-supplied  s^_.cli 
directions  are  desired.  Then  ZANGWL  expects  N  vectors  of  length  N  input  on 
cards  in  4(Fl  5.10,  5X)  Format,  where  N  =  n  is  the  number  of  variables. 

OUTPUT: 

The  output  depends  on  the  print  and  punch  option.  See  tlie  listing. 

ERRORJVIESSAGES: 

None 


PROGRAM  DESCRIPTION: 


ZANGWL  -  acts  as  a  driver  and  convergence  monitor  for  ZGITER.  If 
the  vector  x  =  (X|  ,X2, .  .  .  ,  on  returning  from  ZGITER  is  to  within  ’ 
EACCUR  of  the  value  on  entering,  ZANGWL  returns. 

ZGITER  -  keeps  track  of  the  directions  to  be  searched  and  normalizes 
each  newly  generated  direction. 

POWELL  -  finds  the  minimum  of  the  objective  function  along  a  direc¬ 
tion  supplied  by  ZGITER  using  quadratic  interpolation. 

MATHEMATICAI  METHOD: 

ZANGWL  is  based  on  a  method  proposed  by  W.  1.  Zangwill  in  the 
Computer  Journal,  Vol.  10,  1967,  pp.  293-296.'^  The  method  is  outlined  as 
follows:  Let  f(xj,X2, .  .  .  ,  Xj^)  be  the  function  to  be  minimized,  and  Cj., 
r  =  1 , .  .  . ,  n  be  the  unit  coordinate  directions.  Assume  that  an  initial  point 
p°  and  n  normalized  directions  r  =  1 ,2, . . .  ,  n  are  given. 

To  initialize,  calculate  to  minimize  f(p°+X°t ' ),  then  set  p  ^ ,  = 

’  n  n^n  ’  ^n+1 

Pn  ^n  ^n’  ^  iteration  k  with  k  =  1 . 

k+1  k 

Iteration  k:  p  . } ,  ,  r  =  1, . . . ,  n  and  t  are  given. 

k-1 

Step  (i):  Find  a  to  minimize  f(pf,+j  +  acj).  Update  t  by 


t  = 


t  +  1  if  1  <t  <n 
i  if  t  =  n 


If  a  0,  let  p^  =  pj^^  j  +  crCj.  If  a  =  U,  repeat  step  (i).  Should  step  (i) 


k-1 

be  repeated  n  times  in  succession,  stop;  the  point  p^^^j  is  optimal 


Step  (ii):  For  r  =  1 ,  . .  .  ,  n  calculate  Xj^  to  minimize  f(Pp_|  +  xj^^ 


and  define  =  p|;^_,  +  xj^  •  Let  -  (pj^  -  pj^;‘,  VlIPp  -Pn;) 

k  k  k  vk-I 

Determine  Xj!j^j  to  minimize  f(pj^  +  Xj^^.]  fp+i )  ^nd  set 

_k  -„k,T,k  j.k 

-Pn  +  ^n+l?r 


k  _  .  k  k-1 


-k-1, 


.k+1 


Define  r  -  I , .  .  .  ,  n 


Go  to  iteration  k  with  k+1  replacing  k. 
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NUMERICAL  EXAMPLES  -  GENERAL  TEST  FUNCTIONS 


In  this  section  we  list  the  functions  for  which  ZANGWL  computed  the 
minimum,  and  summarize  the  results.  These  test  functions  are  selected  because 
the  surfaces  they  define  have  steep  curving  valleys  or  have  many  known  optima. 
The  results  are  tabulated  for  each  function.  This  table  lists  the  various  initial 
points,  the  computed  optimum,  the  minimum  value  of  the  objective  function, 
the  number  of  function  evaluations,  and  the  elapsed  CPU  in  seconds’  time  on 
the  IBM  360/65. 

1  2 

In  comparison  with  minimization  routines  (FP  and  CNJGAT)  at  NELC, 
the  computing  times  are  a  bit  slower  and  many  function  evaluations  are  needed, 
but  with  the  speed  of  the  360/65  this  is  relatively  insignificant  compared  with  the 
number  of  man-hours  spent  in  deriving  the  analytic  expression  for  the  gradient. 

1.  Rosenbrock’s^^  function:  f(xj,X2)  =  IOOIxt  -  +  (xj  -  1)^. 

This  function  has  a  steep  valley  along  the  parabola  x-)  =  Xj  with  a  minimum 
at  (1,1). 

2.  Cube:'^  f(xj,X2)  =  100(X2  -  x^)^  +  (xj  -  1 )-.  Cube  is  similar  to 

3 

Rosenbrock’s  function  except  the  steep  valley  follows  the  curve  x-)  =  Xj  and 
has  a  minimum  (1,1). 

3.  Helical;^  f(x  j  ,X2,X3)  =  100  (X3  -  100)^  +  (r  -  1 )“  + x^  where 

n.  ^ 

Xj  =  r  cos27r6,X2=  r  sin27rfl,and  r  =  VX|  +  x^.  This  function  has  a  steep  helical 
valley  with  a  minimum  at  ( 1 ,0,0). 

4.  THREE:^*  f(X| ,X'),X3)  = - ^ - - -sin  ('/27rx-)X3) 

1+(X|-X2)- 

Xl+X2\  \-  ^ _ 

- )•-2 1  .  THREE  has  minima  at  Xi  =  x^  =  xo  =  ±  V  4n+l ,  n>0 

^2  I  I 

integral  with  a  minimum  value  of  -3.  This  function  tends  to  change  quickly 
from  till-  point  (0,1 ,2)  and  then  flattens  out  until  it  reaches  an  optimum.  The 
optimum  depends  on  the  starting  point  Xq. 

5.  FOUR:^“  f(X|,X2.X3,X4)  =  (Xj  +  10x2)“  +  5(X3  -  X4)“  +  (x-,  -  2x3)"^ 
+  10(xj  -  X4)^.  FOUR  has  its  minimum  at  (0,0,0,0). 

6.  CHEBYQ(UAD):^^  This  relatively  new  function  allows  testing  a 
routine  on  a  function  with  an  arbitrary  number  of  variables;  i.e.,  CHEBYQ 


(xj , . .  .  ,  x^),  where  n  is  a  parameter  preset  by  tne  user.  For  n  =  1 ,2, . . . ,  7,9 
the  minimum  value  of  CHEBYQ  is  zero;  however,  for  other  n  the  minimization 
is  still  valid. 

Briefly,  CHEBYQ  does  the  following;  Let  x  =  (Xj, . . . ,  Xj^)  be  a  vector 
(abscissae)  whose  coordinates  are  in  the  range  0  <  Xj  <  1 .  Then,  choosing  the 
shifted  chebyshev  polynomial  Tj,  we  define: 

/•I  1.^ 

^iW=  Jq  Ti(z)dz--2Ti(Xj) 
n 

Then  the  function  fix)  =  ^  (Aj(x))^  has  the  property  that  if  X  is  the  vector 

i=l 

of  abscissae,  then  f  =  0;  otherwise,  f  >  0.  Although  contrived,  CHEBYQ  is  a 
good  example  of  a  complicated  objective  function  that  can  occur.  The 
FORTRAN  IV  listing  of  CHEBYQ  follows. 


SUBRDUTINf:  CHPBYQI  F,X,N) 

IMPLICIT  Rr  AL*«( 

LOGICAL  lEVFN 

DIMENSION  Y(P'')  ,TI(2C),TIMIN(20),X(l) 
DELTA  =  O.ODOO 
ZERO=  O.O'JO'' 

ONE  l.OOOC 

TWO  =  2.0000 

DO  10  J=1,N 

Y(J)  =  rWil*V(,JI-0,NE 

DELTA=  delta  +v(J) 

TK  J)  =  Y(  J| 

10  T^MI^I(J)  =  ONE 

E  =  DELTA  *  DELTA 

IF:VFM  =  .EALSE. 

DO  20  1  =  2,  N 

lEVFN  =  .NOT. [EVEN 

DELTA  =  ZERO 

DO  19  J=1,N 

TTPLOS  =  TKO*Y(  J)<=rnj)-TIM1N{  J) 

delta  =  DELTA  *  TIPLUS 

TIMINIJ)--  TKJ) 

19  TKJ)  =  riPLOS 

A  =  ZERO 

lEIIFVENl  A=  -ONE/( I*I-ONE ) 

DELTA  =  DLLTA/N-A 

20  r  =  F4-  DEL  TA'i'lJFLTA 

RE  TURN 

EF'D 
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ZANCWL-  ROSIE 


optimum  at  (1.0, 1.0) 


QSTEP  =  0.1;  EACCUR  =  1  .OD-04 


Initial  Point 

Computed  Optimum 

Objective  Function 

Number  of 
Evaluations 

Time  (sec) 

1  -1.2 

2  1.0 

O.lOOOOOOODOl 

0.10000000001 

0.16961904 D-22 

325 

0.26 

1  -1.2 

2  -1.0 

O.lOOOOOOODOl 

O.lOOOOOOODOl 

0.16%1905D-22 

328 

0.20 

1  1.2 

2  -1.0 

O.lOOOOOOODOl 

O.lOOOOOOODOl 

0.84351 733  D-24 

206 

0.11 

QSTEP  =  0.01 ;  EACCUR  =  1 .00-04 


1  -1.2 

0.99999991  DOO 

0.7793  5877D- 14 

366 

0.19 

2  1.0 

0.99999992D00 

1  -1.2 

0.99999991  DOO 

0.7793  5878D- 14 

370 

0.19 

2  -1.0 

0.99999992D00 

1  1.2 

O.lOOOOOOODOl 

0.23682822D-14 

152 

0.12 

2  -1.0 

O.lOOOOOOODOl 

ZANGWL-CUBE 


optimum  at  (1.0, 1.0) 


QSTEP  =  0.1;  EACCUR  =  1  .OD-04 


Initial  Point 

Computed  Optimum 

Objective  Function 

Number  of 
Evaluations 

Time  (sec) 

1  -1.2 

2  1.0 

0.99999999000 

0.99999999D00 

0.311345410-21 

402 

0.28 

1  -1.2 

2  -1.0 

0.99999999D00 

0.99999999D00 

0.311348350-21 

399 

0.23 

1  1.2 

2  -1.0 

0.99999999000 

0.99999999000 

0.783690230-21 

179 

i 

0.17 

1 

QSTEP  =  0.01 ;  EACCUR  =  1  .OCMM 


1  -1.2 

2  1.0 

0.10000000001 

0.10000000001 

0.379209870-25 

473 

0.25 

1  -1.2 

2  -1.0 

0.10000000001 

0.10000000001 

0.379209870-25 

470 

0.25 

1  1.2 

2  1.0 

0.10000000001 

0.10000000001 

0.271449010-25 

216 

0.12 

«5 


ZANGWL- THREE 


optima  at  Xj  =  X2  =  X3  =  ±  V'4n+1 ,  n  integral 


QSTEP  =  0.1;  EACCUR  =  1  .OD-04 


Initial  Point 

1 

Computed  Optimum 

(Ajective  Function 

Number  of 
Evaluations 

Time  (sec) 

1  0.0 

O.lOOOOOOODOl 

-0.30000000001 

288 

0.38 

2  1.0 

0.99999999000 

3  2.0 

0.99999999000 

1  0.0 

-0.10000000001 

-0.30000000001 

298 

0.30 

2  -1.0 

-0.10000000001 

3  -2.0 

-0.99999999000 

1  C.O 

0.40012498002 

-0.30000000001 

5710 

4.88 

2  1.0 

0.40012498002 

3  -2.0 

0.40012498002 

STEP  =  0.01 ;  EACCUR  =  1 .00-04 


!  0.0 

2  1.0 

3  2.0 

0.10000000001 

0.10000000001 

0.99999999000 

-0.3000000Ct001 

257 

0.28 

1  0.0 

-0.10000000001 

-O.3C!000000DOl 

257 

0.35 

2  -1.0 

-0.10000000001 

3  -2.0 

-0.99999999000 

2ANGWL  -  HELICAL 


optimum  at  (1 .0, 0.0, 0.0) 
QSTEP  =  0.1;  EACCUR  =  1  .OIMM 


Initifl  Point 

. 

Computed  Optimum 

Objective  Function 

Number  of 
Evaluations 

Time  (sec) 

1  -1.0 

2  0.0 

3  0.0 

O.lOOOOOOODOl 
-0.343  72772D- 15 
-0.81379053I>15 

0.29631 116D-30 

462 

0.56 

1  1.0 

2  1.0 

3  1.0 

O.lOOOOOOODOl 
-0.55742421  D-16 
-0.90526145D-16 

0.77927466 D-32 

366 

0.38 

1  1.5 

2  0.5 

3  A).5 

O.lOOOOOOODOl 
0.86066593D-17 
0.1635661 5D-16 

0.18577895D-33 

351 

0.38 

QSTEP  =  0.01 ;  EACCUR  =  1 .01>04 


1  -1.0 

2  0.0 

3  0.0 

O.lOOOOOOODOl 

-0.54445252D-21 

0.11%3532D-16 

0.14457809D-31 

622 

0.59 

1  1.0 

2  1.0 

3  1.0 

O.lOOOOOOODOl 

-0.50124625D-17 

-0.77793682D-17 

0.6051 8735D-34 

375 

0.42 

1  1.5 

2  0.5 

3  -0.5 

O.lOOOOOOODOl 
0.3741 1888D- 20 
-0.3192021 71>l  8 

0.1 0674 562D-34 

337 

1 

0.33 

ZANCWL- FOUR 


optim’im  at  (0.0, 0.0) 
QSTEP  =  0.1 ;  EACCUR  =  1  .OD-04 


Initial  Point 

Computed  Optimum 

Objective  Function 

Numlor  c*' 
Evaluations 

Time  (sec) 

1  3.0 

2  -1.0 

3  0.0 

4  1.0 

-0.13754371  D-5 
0.13754371  D-6 
-0.5771 5785D-6 
-0.5771 5786 D-6 

0.68469142023 

1350 

0.95 

1  -3.0 

2  -1.0 

3  2.0 

4  1.0 

-0.71682591 1>6 
0.716825900-7 
0.1058991404 
0.10589913  04 

0.36195490018 

1169 

0.73 

1  3.0 

2  1.0 

3  -2.0 

4  -1.0 

0.5815911705 
-0  5815911706 
0.4431292005 
0.4431292105 

0.79920637020 

1568 

0.97 

QSTEP  =  0.01 ;  EACCUR  =  1 .00-04 


1  3.0 

2  -1.0 

3  0.0 

4  1.0 

0.4637473906 

-0.4637474107 

-0.3217002605 

-0.3217002505 

0.35002584020 

781 

0.57 

1  -3.0 

2  -1.0 

3  2.0 

4  1.0 

-0.4071097007 

0.407 10%6O8 
0.4014159306 
0.4014159106 

0.79074373  D  24 

1503 

0.95 

1  3.0 

2  1.0 

3  -2.0 

4  -1.0 

0.2935461606 
-0.2935461607 
-0.45 19863  1 06 
-0.4519863306 

0.36784227  023 

1224 

0.76 

ZANGWL-  CHEBYQUAD 
QSTEP  =  0.01 ;  EACCUR  =  1  .OD-04 


Initial  Point 

Computed  Optimum  | 

Objective  Function 

Number  of 
Evaluations 

Time  (sec) 

n  =  6 

1  0.142 

0.66876591  D-1 

2  0.285 

0.28874067000 

3  0.428 

0.36668229000 

0.320393470-17 

704 

4  0.571 

0.6333 1770D00 

5  0.714 

0.71125932000 

6  0.857 

0.93312341000 

n  =  8 


1  0.111 

0.431 527600-1 

2  0.222 

0.19309084000 

3  0.333 

0.26632870000 

4  0.444 

0.50000000000 

0.351687370-02 

5  0.555 

0.49999999000 

6  0.666 

0.73367129000 

7  0.777 

0.80690916000 

8  0.888 

0.95684724000 

n  =  10 


1  0.090 

0.596199010-1 

2  0.181 

0.16670828000 

3  0.272 

0.23917065000 

4  0.363 

0.39888429000 

1 

5  0.454 

0.39888429000 

0.65039548002 

2603 

14.30 

6  0.545 

0.60111570000 

7  0.636 

0.60111570000 

8  0.727 

0.76082934000 

9  0.818 

0.83329171000 

10  0,909 

0.94038009000 

o  O  r'J  O  O  O  O  O  O  O  O  o  o  r)  O  x*)  O  O  O  O  '■>  o  o 
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c 

r 

C  PUi^TMiSL 

C  TD  FINiJ  THr  '-^KvI^’IM  Qp  A  RFAL  VALUED  FUNCTION  OF  N-VAKMrtLES 

C  WM(!SC  VAL'JES  A^!  E  iJNCHNS  I  » A!  NED. 

r 

r  USAAp 

C  CALL  /A''1C',;L  (  v  I  ,  N ,  E  ACCUK  ♦  OSTEP  » I  STOP ,  LPRI  NT  ,  1  X ,  LPUNCH,  XOPT  ,  FF  ) 

C 

C  nESC>^,  IPT  in.N  OF  PA'^AVFVEKS 

C  X!  -  T,^F  INITIAL  GUESS  OF  THE  OPTIMU'I 

C  'I  -  ■THE  ME  VAPIADLES 

C  EACC'JW-  the  INITT/\i  accuracy  DESIRED.  FOR  REST  RESULTS  SET 

LESS  THAN  OSTEP^^S. 

OSTE?  -  THF  initial  STEP  SIZE  FOR  THE  ONE  DIPENSIGNAL  SEARCH 
POir'TNt  t)OUELL.  USTEP  EQUAL  .1  V^ORKS  WELL, 

ISTCr  -  0STh:P  REfJUCTIGN  CODE 
NO  RE  DUCT  ION 

1  ATTFR  ETwniNG  AN  MINIMUM  TO  WITHIN  FACCURf  QSTEP  IS 
SET  TO  F.ACr.'lR,F  ACCUR=FACC»'R#«2  AND  THE  ROUTINE  DOFS  ONE 
PIML  MIMwizaTIUN  T)  WITHIN  EACCUP. 

I  RRINT-  Pk  tnt  code 

0  no  NOTHING 

1  EACH  N-ni^^PNSlONAL  ITERATION  ONLY 

2  FOMAL  RESULTS  ONLY 
0  1  +  ? 

A  3  OLUS  THE  POINT  ON  ENTERING  AND  LEAVING  POWELL. 
tX  -  CHOICE  OF  I'SER  SUPPLIED  DIRECTIONS.  I  F .  EQ  .C  ,  CO-CRDI  NATE 
OIFrCTIONS  ARE  USED. 

LAUNCH--  A  oiJNCM  COT  ION,  IF  DIFFERENT  FROM  C,  ZANGWL  WILL  PUNCH 
OUTPUT  FOR  the  'nNIMUM  POINT. 

XOPT  -  A  LIST  OP  LENGTH  N  CONTAINING  THE  MINIMUM. 

PP  -  FF  =  FUNC  (XOPT,;N) 

NOTE  ALL  OAPA^ETPPS  ARE  DOUBLE  PRECISION. 

SUBROUTINES  AND  FUNCTlUN  SUBPROGRAMS  REQUIRED 

USER  M,,'ST  S'IPRLY  the  function  SUBPROGRAM  FUNC(X,N> 

METHf'O 

THIS  ROUTINE  IS  BASED  ON  A  METHOD  PROPOSED  BY  W . I . Z ANGW  ILL » 

C  computer  jrinR.,vnL.lOT  1067,  0293-206. 

c 

c  . . . . . . 

c 

SUBPG  JTINL  /ANGWl  {  N^ w  ,  N , E ACCUR  ,  QST £ P  ,  I  STOP,LPRINT,  I X  , l.PUNC H ,  XOPT  , 

1  FF  ) 

implicit  r\E  AL^R  (A-'-,C-M 

DIMRNSIU\  NPkin,  OLDlZuI,  XI(20,2l),  D(20),  PT(20),  XINEW(20» 
DIMtNSlON  R(?C),  XJ(20,2U,  XOPT(l) 

INTEGER  T,  OkIFN,  COUNT 
REAL*R  LAhpDA  ,Ni-  W,  MINEN 

CHMyn-T/z  /'NG/FPSl  LN',Q,E  ,  I  POwEL  ,  OB  JFN ,  L  I  ST,  COUNT 
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IP=0 

Q  =  QST[!P 
C  =  EACCUk 
list  =  LPkINT 
IHMCRD=  1. PUNCH 
ORJFN--0 
IPnwEL=? 

EPSILN=1.0D- 1.5 

C  IP=STOP  FLAG  TERMINATING  SIMULTANEOUS  RFDllCTinN  OF  BOTH  U  AND  E 

C  K^NiJMBER  (JF  CUPRENT  ttfRATIJN 

C  OBJFN  =  CtJRRENT  QUANTITY  OF  OBJECTIVE  FUNCTIilN  EVALUATIONS 

C  OLO( IC)=IC-TM  COMPONENT  OF  ENTER  POINT  FOR  ZGITER 

C  NEW!  IC)  =  ir,-TH  COMPCNENT  OF  INITIAL  POINT  AND  POINT  CCMPUTED  IN 

C  ZGITER 

C  R( !CI=OLC(  IC  )  AT  START  OF  ZGITER 

C  XI(  ICf  !())=  IC-TH  component  OF  ID-TH  NORMALIZED  NONCOORDINATE 

C  DIRECTION 

C  XJ( IC, ID}  =  XI( If  ,  ID!  AT  START  OF  ZGs TFR 

C  XINEW(IC)=IC-TH  COMPONENT  OF  (N+II-ST  ('EXTRA*) 

C  NORMALIZED  NONCOOkO I N A^E  DIRECTION 

C  PT(IC)  =  IC-TH  COMPON'^NT  OF  MINIMUM  POINT  OF  N-0 1  MENS  I  ON  AL 

C  MINIMIZATION 

C  0(!C)-IC-th  COMPONENT  OF  NORMALIZED  NONCOORDI NATE  CIPECTION  OF 

C  ONE  DIMENSIONAL  mj  i ;  aT  I  ON 

C  ALPHA=MINI MUM  STEP  LENGTH  FOR  N-D 1  MENS lONAL  MINIMIZATION 

C  LA;'iBnA  =  MIN  IMOM  PTEP  LENGTH  FOR  ONE  DIMENSIONAL  MINIMIZATION 

C  MINFN=QBJFCT  IVE  FUNCTION  MINIMUM  VALUE 

C  ZGFN=CURkENT  ZGlTFk  ITERATION  OBJECTIVE  FUNCTICN  VALUE 

C  STFN  =  0BJECTIVE  FUNC''‘!ON  VALUE  AT  INITIAL  POINT 

C  J0MP=FIRST  ZGITER  ITERATION  FLAG 

C  nRJ  =  ri6JFN  AT  START  OP  ZGITEP 

C  PTUIF^MATIOUM  NEk<  A^'D  OLD  POINT  COMPONENT  DIFFERENCE 

C 

10  CONTINUE 

IF  (IP.GT.C)  GO  TO  20 
20  IF  (IP.GT.C)  GO  TO 

STFN=FUNC(NEW,^!) 

FP=STFN 
(lRJFN  =  n'3JFN+  1 
GO  Tp  40 

00  STFN=FUMC(PT,N) 

OPJFN=nBJFN+l 

AO  IF  (LIST.  ED.'')  GO  T*'.  qn 

PRINT  2DC,  N,C,r 
PkIMT  43?,  L  I  ST,  IjMCPD,  IX 
IF  (IP.ED.D)  GO  TO  s? 

PKiriT  400,  (  TC,PT(  1C  )  ,  IC-1  ,N) 

GO  Tn  60 

5i;  PRINT  4',^,  (  ir,,\'Pi;(  !C)  ,  IC=l,N) 

60  PR^^If  53D,  SIFN 

IF  (  IP.CQ.''  )  GP  TO  7? 

GO  tt)  so 
7D  Cf;NTINiJE 
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80  IF  (IP.GT.O)  GO  ro  ?30 

IF  (  IX.EO.O)  GO  TO  ^0 
C 

C  USER  SUPPLIES  INITIAL  NHNC OHRO I NA TF  DIRECTIDNS 

C 

READ  340,  ((  XH  IC,  in),lC  =  l,NI  ,ID=1  ,N) 

GO  Tn  130 
C 

C  COMPUTE  INITIAL  NUNC  OGRO  I N  ATE  DIPECTUJNS 

C 

on  DO  120  ID=1»M 

DU  no  IC=1,N 
IF  (  IC.EO.  ID)  GU  TO  ICO 
XI ( ic, ID )=o.onoo 

GU  TO  110 

100  XI ( ir, in)=i.onoo 
no  CONTINUE 
120  CONTINUE 
C 

C  MINIMIZE  IN  ONE  DIMENSION  USING  INITIALIZATION  DATA 

C 

130  00  140  IC=l,N 

140  D(  IC)  =  XI(  IC,N  I 

IF  (IP.EQ.n)  r.U  TO  16C 
IF  (L  IST.EO.r )  GO  TO  145 
PRINT  580,  (  tC,PT(  ir  niC=l.N) 

145  CALL  POWELL  (  N  ,  OLD  ,  I) ,  LAM60A ,  FF  ) 

DO  150  IC=1,N 

150  OLOI IC)=QLO( IC ) fLAMRnA*0( IC) 

IF  (LIST.EO.O)  GO  TO  155 
PRINT  8^0,  nC,nLO(  IC),IC=1,N) 

155  CONTINUE 
GO  Tn  18C 

160  IF  (LIST.lQ.O)  go  in  lo5 

PRINT  560,  nC,NEW(  IC  »  ,  IC=1,N) 

165  CALL  POWELL  ( M , N EW , 0, LA mboa  ,FF  ) 

C 

C  COMPUTE  FIRST  POINT 

C 

DO  no  IC=1,N 

170  ULI)(  !C)=Nrw(  I  C  )  ,-i  4MH0A*n(  IC) 

IF  (LIST. to. n)  GU  TO  175 
PRINT  ^00,  (  ICjOLPI  ’C),1C=1  ,N) 

175  CONTINUE 

r, 

C  INITIALIZE  FOP  E  IP  ST  nf-VATIUI 

C 

K=  I 
T=1 

JUMP= 1 

180  IF  (ix.ro.^l  GO  Tf;  m 
JUMP- JU'MPf  I 
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19'  DO  ’’(C  IJ=1,N 

200  k(  f  j)=r)Ln(  I  j) 

nr  720  iK=i,N 
DO  210  IL=1,N 
210  XJ(  IL , IK  )  =  X1  (  II  , TK) 

2?r  COMTTmijo 
npj  =  nHjFM 

IF  (LIST.NF'.M  on  T'l  ^^0 
I'PINT  540,  K 
C 

230  CALI  20ITLL'  (  9  ,0  LU  ,  X  I  ,  r ,  AL  PHA  ,  f  T  ,  X I  NE  W  ,  J  OMO ,  nC)  YE  S  ,  FF  ) 
r 

?/,C  Zr,Fn  =  FF 

c 

C  COMPARE  Nl.F  Ar;n  01. n  ‘MIINTS  FOR  ,'1-0 1  MF  N  S I  OMAL  MINMOM 

ACHIEVED 
C 

PTnrF  =  nA!5S(PT(  1  )-nLO(  I)  ) 
on  2  50  IJ  =  2,,\' 

()M  =  f)AiIS  (PT(  IJ  l-riLIX  IJ)  I 
IF  (ON.LE.PTnrp )  on  TO  2 5C 
PTDI F=CM 
2  50  CGOTINiJF 

IF  { i'>TlHF  .(iL-  .r  )  r.i)  rq  2hC 

IF  IK. lit.  2)  GO  ro  200 
UBJFM-rOHJ 
K  =  K-l 

2  60  IF  (LIST.FO.J)  go  Tf)  7  7C 
IF  ILIST.Fn.2)  GO  TO  270 
PRINT  430,  (  IC  ,R  I  IC  )  , IC=1 ,N) 

PRINT  450,  I  (  IC,  IO,XJI  IC,!())  ,IC  =  1  ,N)  ,Il)  =  l,N) 

PRINT  550,  K 

PRINT  431 ,  (  K  ,ni.i)(  IC  ) ,  IC-  1  ,N) 

PRINT  450,  ((  TC,  ID, XMIC-l!))  ,10=1  ,N)  ,11)^1, N) 

PRINT  510,  I  IC,XINCN(  IC? ,IC  =  1,N) 

PRINT  520,  K,.7GFN 

PRINT  46C,  ( IC ,P  T(  10  )  , IC  =  l  ,N) 

PRINT  ^20,  K,nRJFN 
IF  (PTniF.LT.T)  GO  TO  ZS: 

0 

G  TEST  Fil-  ''-niH^NSltilAL  MLNIMiJv  AOHIFVFD 

2 T'l  IF  (OnUNT.LT.N)  0,1  TO  320 
C  N-D  TMCMS  in.NAL  ACHiEVED 

0 

TRO  Cn.NTINUF 

IF  (  I3MCK0.LE  .  1  )  Gi!  T^’  2'?^^ 

PHNC'^  MC,  I  >,  I  X 
PUNCH  37C,  -i.OjE 
PUNCH  34  0,  ( PTI  IC  !  , 10  =  1 ,N) 

PUNCH  410,  L  1  ST,  N?MrR!),KI,!)CK 
If  IIX.E'J.T)  GO  TO  2'0 

PUNCH  3,C,  (  (  X  I  (  IC  , ’0  ) ,  IC=  I  ,'N)  ,  10=1  ,N) 

PG"  CUNTIN JF 
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M!NPM=  TF 

IF  (L  isr.L  F.  1 )  ^.ri  in  ^ic 
PFIMT  350,  MINFM 
PF-  I  NT  ^60,  (  ir.PTI  TO  ,If  =1,MJ 
PRINT  ^40,  r!0jFN,R 

r.  TFST  PiJK  FINISH  ANN  RESTART  IF  NOT 

310  IF  (  IP  .r.E.  I  STOP  )  GU  TO  3  3C 
Q=F 
F  =  E*r. 

K=K+  1 

JliMP=JIJ‘1P+  I 
1P=  IP+1 

IPnwFL=lPC)WFL*IPni^tL 
GO  TO  10 

C  REITFPATt  U'jrTF  PEG'JlPrf)  ACCURACv  ACHIEVED 

3?o  K=K+1 

JlJMP  =  JilMP  +  1 
GO  TO  IQC 

330  DO  33S  J=1,M 
335  XOPTI J )=OLn( J ) 

PRINT  6?C 
RETURN 
C 

340  FORMAT  (4(1-15.1.0,5X11 

350  FORMAT  (  HO  ,  3  3H03  JFf  T  I  VE  f  UNCTION  MiNIMlJM  VALUE  =  ,  E  13 . 1 1  » 

360  FORMAT  (lHr,2lHM0.  MINIMUM  PDI N T/ ( I  3 ,3 X , El  8 .  1 1  )  ) 

37C  FCjRMAT  (  I<l,?X,ri4,  io,3X,F15.  ICI 

390  FORMAT  (lUl.lfiHNC.  CF  DIM  OR  VAR= , I  2 ♦ 5 X , 16  HONE  DIM  START 
(J  =  ,E18.1l  1  ,5X,9MACCURACY=,E18,  III 
4C0  FORMAT  (lM0,2lHNn.  INITIAL  Pf)  1  NT/ (  1 3  ,3  X  ,  E 1 8.  1 1  )  > 

410  FORMAT! 2 15) 

420  FORMAT  (  1 MC  ,  2  1  m  I  TFP  A  T  T  UNS  I  THROUGH  ,19,10(1  REQUIRED  , 

I9,31H  OflJFX  HIVE  FUNCTION  EVALUATIONS) 

430  FORMAT(28ll  ENT[P  /GITF.R  WITH  THE  PUINT,/4U  NO.  ,  /  ( I  3 , 3X  ,  E  18 . 1 1  )  > 

431  FI1PMAT(26H  LIAVE  /GITCR  AT  THE  POINT, /4H  NO.  /  (  I  3  ,3  X  ,  F  1 8 , 1 1  )  ) 

440  FPR’-'AT  (  IHC,  32HACHIEVING  THIS  MINIMUM  REQUIRED  ,19,361) 

OBJECTIVE  f  IUNCTIUN  EVAUjATTHNS  AND  ,I9,11H  ITERATIONS) 

450  FORMAT  (1H0,31H  NU.  XI  NO,  D I RFO T I  UN/ ( 2X  ,  1 3 , 5X , 

! 3,4X,E18  1.  1  1  n 

460  format  (1HG,17HNC}.  NEW  POI  NT  /  (  I  3 , 3  X,  E 1 8 . 1 1  )  ) 

48C  format  (  IHO,  1  1  wpr  imt  r.UDE=  ,  I  3 , 5X  ,  1 1  UPUNCH  CODE  =  ,  1 3 , 5  X  ,  1 2HT  I  M  I  NG 

CU1UE=,  !3,5X,PHXI  f  LAG=,I3» 

51C  format  (lHr,?4UN0.  XI(N+1)  / ( I  3 , 3X  ,  E  18 . 1  1 )  ) 

5?n  format  ' IMOv IPHT TFRATIUN  ,I9,37U  CHANGES  OBJECTIVE  FUNCTION  VALUE 
1  TO  , PIR.  1  1  ) 

5?D  FORMAT  (  IHC  ,  ARHOBjEC  T  I  VE  EUNCTIOM  VALI.IF  AT  INITIAL  PP  I  NT  = ,  F  1  8 .  1 1  ) 
540  format  (  lH0,?3iiFNTtF  /GITFP  ITERATION  ,19) 

550  FORMAT  (  1 1 10 , 2  3  HI  FA  VF  ZGITfR  ITERATION  ,1Q) 

5t>0  FORMAT  (lHr,32MNn.  f  M  I  T  I AL  I  Z  A  t  l  f)N  ENTER  POWEl  L/ (  I  3 , 3X  ,  E  18 . 1 1  )  ) 
580  FORMAT  (1U^,17M^().  CNTFR  POWELL  /  ( I  3 , 3  X  ,  R 1  8 .  1 1  )  ) 

5Hu  format  (1U!M7hNQ.  leave  PnwCLL/(  I  3,3X,E18. 11  )  ) 

60D  FERMAT  (  1  ,  IRHND .  INITIALIZATION  LEAVE  POWELL/ (  I  3 , 3x  ,  F  18 . 1 1 )  ) 

610  FC'RMAT  (211P) 

62^  FORMAT  ( lUl ) 

EM  D 
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SUBPnUTINE  ZGITFR  (  NOV  i HPT  ,1)  I  R  ,  NC  f  DI  ST  ,  G ,  P X TP  A  ,  J  SW  ,  I  SW  tPP) 

IMPLICIT  PFAL*fl  (A-H, 0-7.1 

CIMFNSIllN  nPT(?G),  niR(20,2n.  C(20)»  G(20),  0(201,  r!rw(2'^),  LXTRA 
I  (20 

RFAL*R  L,NCW 
IMTFGEP  niUTN,  COUNT 

CnMM0N/7ANG/t-PS  1LN,0,F  ,  I  PUWCL  ,()OJFN  ,L  I  ST, COUNT 
r 

C  UNCONSTRAINLD  M QV- U I  MENS lONAL  wimtmi/ATION  WITHOUT  USING  UFPIVATIV 

C  USING  GIVEN  POIN^  DPT  AND  GIVFN  DIRFCTKJN  OIK 

C  NPV  =  OUANTITY  OF  VARIARLFS  IN  OliJFCTIVfc  FUNCTION 

C  MC=MUMRER  Of-  CURRENT  COORDINATF  PIRECTION 

C  DPT(  IC)=  IC-TH  CriMPHNENT  OF  OLD  POINT 

C  NEvM  IC  )=  IC-TH  COMPONl^NT  OF  NEW  POINT 

C  DIR  (  IC  ,  ID  )  =  IC-Ti(  COMPONENT  OF  !  D-TO  NONCOOKO  IN  AT  E  ki1RMALI7t:o  DI«EC 

C  FXTRAI  I  C  )=  IC-TM  COMPONFNT  OF  (NUI-ST  ('EXTRA') 

C  M0RMALI7h;D  NOMC  OOP  0  I N  AT  F  DIRECTION 

C  OlST=MINIHilM  STE’’  LENGTH  ALONG  COORDINATE  DIRECTION  NC 

C  C(ir)=IC-TF)  CO^'PONFNT  of  CURRENT  NURMALIZED  COriRDlNATF.  DIRECTION 

C  G(IC)=IC,-TH  component  of  MINIMUM  POINT  IN  NDV-MI  M  MI  ZAT 

C  L=MINIMUM  STEP  LENGTH  IN  UNL  DIMENSIONAL  MINIMIZATION  ALONG 

C  NONCnOROINATE  DIRECTION  H 

C  H(lC)  =  IC-TH  COMPONENT  UF  CURU^ENT  NONCUOPf)»  NATe  DIRECTION  IN  ONE 

C  DIMENSIONAL  MINIMIZATION 

C  CnilMT  =  CURr<  Ent  TOTAL  CE  COOPDINATE  DIRECTIONS  USED  WItrOOT  COMPOT  IN 

C  'EXTRA'  Ni'iNCOORO  ■^lATE  DIRECTION 

C  JSW=riRST  ZGITcp  ITEPATIOM  FLAG 

C 

C  PART  ONE 

COUNT  =  0 

C  IF  K^l  AND  IX  (IN  ZANGWILL)  =C ,  GO  TO  '>ART  TWO 

IF  (JSW.GT.l)  GO  TO  2G 
NC=NC+1 

DQ  lO  JK=1,NDV 
ID  MEW( JK )=OPT( JK  ) 

GO  Tfi  UG 
C 

C  COMPOTE  CU'^RENT  ClilPn  TnATE  DIRECTION 

C 

20  no  ^,0  JC  =  1,NDV 

IF  (JC.EQ.NC)  GO  rn  00 
C( JC)=D.ODOO 
GO  TO  AO 
c(  jc)  =  i.':d().^ 

40  CONTINUE 

C 

C  MINIMIZE  IN  i1NF  DIMpnSIIN  ALONG,  CURRENT  CfHDROINATF  niRFCTI';\ 

C 

IF  (LlST.NE.i.)  Gn  TO 

PRINT  2S0,  ( JC ,DPT( JC  I  , JC=  1  ,NnV) 

4S  CALL  PUWtLL  (  NT)  V  ,0  P  T  ,  C  ,  D  I  S  T  ,  F  F  ) 

C  TLST  FOR  ALL  COORD  IMA'^E  DIPECTIONS  UGtO 

IF  (NC.NE.NOVI  '',0  ffi  OC 
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N'C=  1 

no  TH  60 

sc  NC=NC+1 

C  TiST  FDK  MINIMUM  STi^P  LtNGTH  ALUNG  C(NC) 

SO  If  (TABS  (  MIST)  .r.L.FOSILNI  GO  Tf)  OQ 

1  r  (  Cl  ION  T  .  or. .  NO  V  )  GI  i  to  70 
Cl.UMT  =  Ci)UNTf  1 

on  Tf)  20 

C  MII'IIMDM  STCP  LTNCTH  ALONG  f.  (  NC  )  f-Gi.IND 

70  DO  HO  ,JK^I,N[)V 

SO  G(  JK  )  =  n'->TIJK  ) 

fftmhn 

c 

C  UOIJATL  CO'^HfNT  I’MINT 

90  DO  KO  JC=1,NI'>V 

lOr  NCWI  JC)  =  Df»T(  jf  )  tniST^CI  JC) 

IF  (L  IST.nF  .S  )  GO  Tr)  1  IC 
(-'F’NT  2(.^r,  (  jr,!JFv,'{  JO  ,  JOl  .NDV) 

r 

C  PAPT  TWM 

0 

C  KINI'^IZE  IN  HNf  DIMENSION  USING  CORPFNT  POINT  AND  CDPREN 

C  NONC!iUP,niN/.TC  iJliTXTIONS 

lie  no  iHc  jn=i,Nnv 
DO  120  jc  =  i,r4nv 
120  HOC  )  =  ')IW  (  JC,  JO) 

c 

C  TEST  F'lK  7I;fC!  OIRFC.tion 

I7=C 

no  no  j=L,N':iv 

IF  (t)A  '.S  (H(  J)  )  .LC. 1.00-15)  GO  TO  130 
17  =  1 

no  contimuc 

IF  (TZ.HC.O)  GO  T't  in': 

IF(L  IST.r'L  .A)  GO  TO  155 

PRINT  250,  (  IJ,NPW(  IJ),lJ=l,NnV) 

155  CAlL  f'OWFLl  INOV  ,NFW  ,H,L,FF  ) 

ibC  on  no  K  j  =  I  ,Nnv 

1  TO  NF  (  KJ  )=Nt;W(  K  J  I  K  J) 

IF  (1.  IST.NF  .4)  GO  Tfi  nO 
PKI^'T  2^,^,  (  J  1  ,NF'N(  J  I  ) ,  J  1=  1  ,NDV) 
lOG  contindf 

C 

C  COyP'ITE  'PXTt.'A'  NONCnORD  INATE  DnECTlCN 

C 

DFMOM^O.  ODG'^ 

Do  no  Jf,=  IjNDV 

inr^  Of  NOM^  JLNOMfOAHS  IN  F  W  (  JC  ) -1)p  T  (  JC  )  ) 

DO  200  JC=1,NPV 

DIR  (  Jf  ,NDV*  I  )  =  (  fD  W(  JC  »-JPT  UCI  )  /nF:NQM 
?^0  E  XTP  A(  JC  )  =  ')T‘-'  (  JO  ,Ni)V>  1 ) 
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c  MINIMIZE  IN  0^‘r  niMENSIdN  AL()^4G  'EXTRA*  N'lNC'inKD  IN  AT  F  DIRFCTION 

1)'’  ?10  JC=1,MDV 
210  HI  jr.)=OIK  (  JC  ,Nnv+l » 

IF  (LIST.NF.4)  ON  TO  ?15 

PRINT  2<so,  (  JC  fNF.WlJC  )  f  jr.  =  l  ,MI)V» 

?15  CALI.  PGWI  i.L  IN[.iV,NFW,HtL,FF  I 

nn  jc=i,Nnv 

2  2^  DPT  IJC  )  =  NFN(jr  UL*H(  JC  ) 

IP  (1.  IST.NF.4)  GO  TM  ?2‘) 

PRINT  ?60,  (JC,NPT(  jr.  l,JC=l,Nl)V) 

C  CfiMPl.lTF  MHNCNf'P. D  1 N ATT  DIRECTIONS  FOR  NEXT  ITf.RATIQN 

22S  1)0  2^^ 

DO  2'T  JC-l,NDV 

DIP!  JC,JO)=DFI  JC,JD+1» 


CONMNOr 

140  DO  ISO  JK=1,NPV 
Gl  JK  )  =  NF'aI(  JO 
P  FTORN 
r 

7S0  FORMAT  I  up:  ,  1  thNO.  i'ORTLI  E  NT  L  R/ ( I  .S  ,  3  X  ♦  C  1  B  .  1 1  )  ) 

240  FORMAT  (IHO.ITHNO,  POWELL  L E A VE / ( 13 , 3 X , r 1 s . O ) ) 

END 


S'lHP/lLITINF  POl-MLL  (  N  , '’0  •  S  t  S  TL  P  t  F' F  ) 

U-‘PI  'CIT  RPAL=!-'«  I  A-'Un-Z  » 

DIMENSION  SI?''),  D0I?0),  ?{?Q),  F{4),  VI4) 

INTEOFP  OfijPN,  COUNT 

COMMnN/ZANO/FF  SILN  ,0,F  ,  1  PNWt L  , O')  JFN ,  L  I  S T ,  C HUNT 
C  UNCONSTF  AINFi)  ONP  N  I  MP.NS  IONAL  M  I  N  I  M  I  Z  A  )  I  <')')  WlTHfUJT  OST'NO  DERIVATIV 

finding  OFLOW  reasonable  Q  value,  see  ;  CALAHAN,D.  A.  ,  COMPUTER 
C  A!i)ED  network  design.,  MCGRAW-HILL,  1948, 

C  P0=  INITIAL  POir:T  VECTOR 

C  S=G!VFN  DlkECTION  VICTOR 

C  GIVEN  Pn,S;  FINN  STEP  THAT  MINP-U/tS  OBJECTIVE  FUNCTION  PO+STEP*S 

C  GIVEF)  PO,S.  flNO  STEP  THAT  MINIMIZES  OHJECTIVP  F  UNCT  lONI  pC*STE  P* 

C  N^UOANTITV  of  VAwIARLFS  in  objective  FUN^TI^)^J 

C  (U  INITIAL  STEP  LENGTH  Al  llNG  S 

F'  F  =  GEOMlTPir.  S^P'PS  p'ATIG  PDF'  FINDING  REAS'JNABLE  STPp  LENGTH  IN 

C  G'TMPUTUiG  SF  T  OF  INITIAL  THREE  PUINTS 

0  CF  =  Gt;nMETR  IC  ^FRIFS  COEPFICIENT  FUR  FINDING  ABOVE  RFAS^^,"BLE  STFP 

C  E=REOllIP,ro  ACCURACv  UF  MJNIMijm  POINT  (FACH  COMPONENT) 

C  V=ARKAY  OF  CUPPFr.'T  PfJINT  VALUES 

C  F=APPAY  'Or  rjlPRENT  pO'N^  fjBJECTIVE  FUNCTION  VALUES 

C  OBJ  FN=0UANT  I  TY  nr  OPJ.^TiVt  FUNCTION  EVALUATIONS 

G  NUM=Oi|ANT!  ty  'IF  jOADRATIC  I  NTE  R  POL  A  T I  UN  S 

NOM-n 

Jj  COMPuiL  T'lPPl  STARTING  VALULS  AND  THEIR  OBjrr.FlVE  FUNCT'Qri  VALUES 

1  .  40P0 
A=O.PDOO 
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on  1=1 tM 
10  z(n  =  PO(n 

FA^f-P 

UO  20  1=1, N 

20  Z  (  I  >=^P0(  1  I  +0’4‘S(  I  ) 
rC=FONC( Z,Nl 
UHJFM=OnjFN+l 
IF  (FQ.LT.FAI  QJ  Til  *^0 
DO  30  J  =  1,N 

30  Z( I )=P0(  I )-0*S(  1) 

FMQ=PUNC ( Z . N  ) 
UHJFN=0njFN4-l 
IF  (FNQ.GE.FA1  GU  TO  40 
G=-0 
hR=FNO 
ISW  =  n 

GO  TO  60 
40  V(l)=-0 

V( 2  1=0.0000 
V( 31=0 
F(  1  )=FNQ 
F( 21=FA 
F( 3)=FQ 
GO  TO  120 
0  0  0=0 

FP=FO 

ISW=l 

60  CF=i.cnno 

SU^=  l.onoo 

70  CF=CFX‘R 

suM=suM+cr 


IF  (  ISW.EO.  I  I  GO  TO  eO 
C=-0*S0M 
GO  I'n  00 
so  C=0*SUM 

90  on  100  I=1,''M 

100  Z(  I  )  =  P0(  I  )  I  1 

FC  =  FilMC  (  Z  ,M  1 
C!RJFN  =  ()njFN+l 
IF  (FC.OT.FOI  ou  in  lie 
A=B 
FA=Pil 
U  =  C 
FB=FC 
GP  TO  70 
110  v(  n  =  A 

V(21=R 
V( 31=C 
F(  ri  =  FA 
F( 2  J  =  FB 
F(  3)  =  pr. 

c 
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C  COMPUTE  ;:hRn  of  first  OFRiVATIVf  OF  AR  PKOX I  MAT  I  NO  QUADRATIC  THROUG 
f.  THREE  CJIRRENT  POINTS  AND  ITS  OBJECTIVE  FUNCTION  VALUE 

C 

1?0  W1  =  V(2)-V(3) 

K2=V(2)+V(3I 
W3=V(3)-V( 1 » 

W4=VI 3)+V( 1 > 
kS=V(  n-V(  21 
W6=vn  l+V<  2) 

W7  =  2.<'(W1*EI  1  )*W3*F(  2H-WS«F(3)  I 
V(4)=IW1<'W2>!'F(  1  )+vJ3*W4«F{2»*-W5*Wft*F(3I  ) /W7 
Dfl  130  1  =  1,  N 

13C  Z( 1 )=P0( I )+V( A)*S(  n 
r(4)=FUNC(7,N) 
nRJFN=OBjrM+l 

IF  ( V(4! .NT  .V( 7  n  GO  TO  140 
GO  TO  ICC 

C  TEST  THREE  CUPPENT  POINTS  FOR  CLOSENESS  TO  ZERO  OF  FIRST  OERIVATIV 

140  IJ1=1 

ISO  Yl=nAHS  (V(IJ1)~V(4)I 

IF  (VI- e.GE.C.odoC  )  0(J  TO  170 

C  ZERO  OF  FIRST  OFPIVATIVE  IS  MINIMUM  DISTANCE  ALONG  S 

140  STFP=V(4) 
f F=F(4) 

IF  IF(4)  .LT.riOn  GO  TO  240 
STFP  =  V(2  ) 

FF  =  F( 2  ) 

GO  TO  240 

170  IF  njl-2)  180, ISO, ■’00 
180  IJl=2 

GO  TO  15C 
100  IJI=3 

GO  TO  150 
C 

SHRINK  CLOSED  INTERVAL  CONTAINING  mimmuM  BY  DISCAPDING  BOTH  POINT 
nUTSIOF  OF  interval  and  its  OBJECTIVE  FUNCTION  VAlUE  ANC 
RELABELLING  REMAINING  POINTS  AND  THEIR  OBJECTIVE  FUNCTION  VALUE 

700  IF  (V(4)-V(2I  .GT.c.rnoCI  0-0  TO  220 
ir  (M4)-F(  2)  .GT.O.CODGI  GO  TO  210 

F(31=F(2) 

V(3)=V(2) 

FI2  )=F(4) 

V(2)=V( 4) 

IE  (NUM.GE  .IPOWFL  )  F-O  tq  i  ec 
NUM  =  NUM+  1 
GO  TO  120 
210  F(1)=F(4» 

Vn  )=V(4) 

IF  (NUM.GF  .  I  POWFL  I  TO  160 
NUM  =  NUM+  1 
GO  Tn  120 
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IF  ( F(4)-H( ?» .1  r.O.  »  GO  TO  220 
l-(3)^F{4) 

V( 3)^V<4) 

IF  (  NUM.r.b'  .  IFI'WCL  )  GO  Fl)  160 
Nl)F'  =  NlJM4-l 
GO  T(i  12C 
2  20  F ( 1 )  =  F ( 2 ) 

V(  n=V(2) 

F(?)=F(4) 

V(  ?)='/(  4) 

IF  (  NUM.Gt. .  IPnwFL)  on  TO  160 
,MUM  =  NUM^  1 

on  TO  120 
240  CONTIMUF 
FFTIJi^N 
FMO 
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NELC  FESDIR 


CATALOG  IDENTIFICATION: 

E4  FESDIR 

PROGRAMMER: 

Gail  Grotke,  Decision  and  Control  Technology  Division 
PURPOSE: 

To  minimize  a  function  f(Xj,X2, .  . . ,  Xj^)  of  n  variables  whose  values 
are  constrained 

RESTRICTIONS  AND  LIMITATIONS: 

The  function  and  constraints,  and  their  first  partial  derivatives,  must  be 
continuous. 

LANGUAGE: 

FORTRAN  IV 

COMPUTER  CONFIGURATIONS: 

IBM  360/65 

ENTRY  POINTS: 

FESDIR 

SUBPROGRAMS  AND  WHERE  REFERENCED: 

Programmer-Supplied  Programs 

FUNC  called  by  FESDIR,  and  POWEL,  and  BNDY 
POWEL  called  by  FESDIR 
BNDY  called  by  POWEL 

Library  Subprograms 

SORT  called  by  FESDIR 

ABS  called  by  FESDIR,  POWEL,  and  BNDY 

DEFINITION  OF  VARIABLES 


xo 

Initial  feasible  point 

X 

Working  point  and  final  minimum 

s 

Direction  vector 

GF 

Gradient  of  the  function 

C 

The  constraints 

GC 

Gradients  of  the  constraints 

DELG 

Normalized  gradient  of  the  violated  constraint 

DELF 

Normalized  gradient  of  the  function 

NUM 

Number  of  iterations  through  Powel 

N 

Number  of  variables 

IC 

Number  of  constraints 

STEP 

Value  returned  from  Powel  that  gives  the  minimum  in 
the  S  direction 

IFLAG 

Number  of  violated  constraints 

FBOIJND  Criterion  by  which  it  is  determined  whether  function 
values  are  converging 

GBOUND  Criterion  by  which  it  is  determined  whether  direction 
vector  is  converging 

CODE  Code  that  determines  what  will  be  printed 

INPUT  FORMAT: 

A  driver  and  a  function  subprogram  are  needed. 

The  driver  calls  FESDIR  (F.  XO,  N,  FBCUND,  GBOUND). 

The  parameters  N,  IC,  FBOUND,  GBOUND,  CODE  and  the  feasible 
point  XO  lor  FES  DIR  ana  the  parame'ers  E,  IPOWEL,  and  OBJFN 
for  POWEL  must  be  set  in  the  driver.  These  three  parameters  for 
POWEL  are  in  the  common  block  labeled  ZANG.  IC  and  CODE  are 
in  common  with  FFSDIR. 

The  function  FUNC  has  a  parameter  list:  (K,  X,  C,  GF,  GC). 

If  k  =  1  when  FUNC  is  called,  only  the  value  of  the  objective  function 
is  returned. 

if  k  =  2,  the  cons*raints  are  evaluated  and  returned  in  C. 

If  k  =  3,  the  gradients  of  the  constraiiit  and  the  gradient  of  the 
function  are  returned  in  GC  and  GF,  respectively. 

If  k  =  4,  only  the  gradient  of  the  function  is  returned. 

OUTPUT  FORMAT 

Prints  out  according  to  the  Print  Code  CODE 


ERRORMESS.<*.G^: 

None 

PROGRAM  DESCRIPTION: 

Main  Program 

The  driver  sets  the  values  of  the  parameters  and  calls  FESDIR. 

Subroutines  and  Functions 

FESDIR  -  finds  the  minimum  of  the  function  within  the  constraint 
set. 

POWEL  is  a  one-dinicnsional  quadratic  search. 

FUNC  evaluates  the  objective  function,  the  constraints,  and  the 
gradients  of  each. 

BNDY  checks  to  see  if  a  point  chosen  by  the  quadratic  search  is 
within  the  constraint  set  ind  records  the  number  of  con¬ 
straints  that  are  violated. 
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3|2|1  0|  Value  of  CODE _ 

Function  value,  ooint,  and  constraint  at  each  step  in  POWEL _ 

\/\  y/  Function  value  returned  from  POWEL,  step  size,  and  minimum  point 

v/'v/l  Gradient  and  normalized  gradients  of  violated  constraints  and 

function,  and  their  sum  _ 

\/  \J  \J  sj  Initial  and  final  points  minimum,  number  of  iterations  througli 
POWEL,  and  number  of  function  evaluations 


MATHEMATICAL  METHOD: 

Given  an  objective  function  f(x  j  ,X2, . . . ,  x^)  and  a  iet  of  constraints 
gifX]  ,X2, .  .  . ,  X|^)  <  0,  FESDIR  finds  x  such  that  f(x)  is  minimum  and  each 
constraint  gj(x)  is  satisfied ;  that  is,  less  than  zero.  To  accomplish  this,  FESDIR 
uses  a  method  of  feasible  directions  given  in  the  article  “Nonlinear  Programming 
with  the  Aid  of  a  Multiple-Gradient  Summation  Technique”  by  Klingman  and 
Himmelblau.^^  Klingman  and  Himmeiblau  suggest  using  a  new  direction  given 
by  (NSD),  new  successful  direction,  defined  as 

V  Grad  Fix)  I 

(NSD)  =  >  < - - -  + - , 

(iGradCj(x)l  |Grad  F(x)| ) 

where  KC  is  the  number  of  constraints  violated  and  Grad  is  the  gradient  or 
first  partial  derivative.  So  (NSD)  is  the  sum  of  the  normalized  gradients  of  the 
violated  constraints  and  the  normalized  gradient  of  the  function. 

FESDIR  uses  this  idea  in  choosing  a  new  direction.  Starting  with  a  feasible 
point  -  that  is,  a  point,  X,  that  satisfies  all  the  constraints  -  and  with  a  direction, 
S,  FESDIR  uses  a  one-dimensional  quadratic  search,  POWEL,  to  find  a  value 
STEP  such  that  the  minimum  feasible  point  along  the  direction  S  is  given  by 
Xmin  =  X  +  STEP  *  S.  If  any  constraints  were  violated  in  finding  the  minimum, 
a  new  direction  is  determined  by  the  negative  of  the  sum  of  the  normalized 
gradients  of  violated  constraints  and  the  normalized  gradient  of  the  function. 
Otherwise,  the  direction  is  the  negative  of  the  normalized  gradient  of  function. 

The  point  and  direction  are  then  used  to  find  a  new  minimum  point  and  a  new 
direction.  This  process  is  continued  until  it  satisfies  the  convergence  criterion; 
that  is,  until  the  function  values  converge  or  the  direction  vector  converges. 

1 .  The  Two-Variable  Problem 

Minimize  f(X|  ,X2)  =  (xj-2)“  +  (x-)-!)^ 
subject  to  gi  =  x  j  +  X-,  -  2  <  0 

§2  ~  ^  1  ~  x-j  ^  0 

Minimum  f(l  ,1 )  =  1 

2.  Tlie  Circle  Problem 


Minimize 

((Xj.Xt)  - - 

(Xj-  1 )-  +  Xt 

Subject  to 

g|  =-x^-x;  +  4<0 

■>  1 

g-)  =  -16  +  x“  +  XT  <  0 

Minimum 

1 

O 

II 

b 
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3.  The  Three-Variable  Problem 


3  2 

Minimize  f(x  j ,  X2,  x^)  =  x  j  -  6x  j  +  1 1  x  j  +  x^ 

2  2  2 

subject  to  g|  =  Xj  +  X2  -  Xg  0 

_  2  2  2^... 

g2  =  -Xj  -X2-X2-(-4<0 


g3  =  -x,  <0 

84  =  -X2  ^  0 

g5  =  -X3  <  0 

g6  =  X3  -  5  <  0 

Minimum  f(0,  V2~,  V^)  = -V^ 

4.  The  Colville  Problem^^ 

1 .  The  Two-Variable  Problem 


Initial  Point  Computed  Optimum  Objective  Function 


1.  -1.0  0.99999982  1.0000003 

2.  2.0  0.99999969 

2.  The  Circle  Problem 

1.  2.0  -0.199589  EOl  -0.988216 

2.  2.8  0.141888  EOO 

3.  The  Three- Variable  Problem 

1.  0.378964  EOO  0. 1 7968266  E-08  1.5673426 

2.  0.168076  EOl  0.15673425  EOl 

3.  0.234720  EOl 


4,  The  Colville  Problem 


1.  0.78619999  E02 

2.  0.33439999  E02 

3.  0.31070000  E02 

4.  0.44180000  E02 

5.  0.35319999  E02 


0.78012675  E02 
0.33144035  E02 
0.30144462  E02 
0.44999996  E02 
0,36546304  E02 


-0.30631831  EOS 


Number  of 
Iterations 

32 


191 


30 


29 
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SAMPLE  PROBLEM  AND  PROGRAM  LISTING 


C  FEASIBLE  DIRECTIONS 

DIMENSION  XCS),  XO(5» 

INTEGER  OBJFN.CODE 
COMMON/ZANG/  E ♦ TPOWELt OBJFN 
COMMON/CON/ICf KQ, CODE t ICON ( 16) 

KQ=6 

CODE=l 

N=2 

IC=2 

XO(  1)=-1.0 

XO( 2)=2.0 

FBQDMD  =  l.OE-12 

GBOUND  =  l.OE-12 

E=1.0E-6 

IPOWEL=5 

OBJFN^O 

CALL  FESOIR  » F , XO , X , N, FBOUND* GBOUND » 

STOP 

END 

FUNCTION  FUNC  ( K , X ,C , GF , GC ) 

DIMENSION  C(  16)  ,GF(5)  .GCn  6,5)  ,X(  5) 
FUNC  =  1.0 
TA=X(  1  )-2.0 
TB^X(2)-1.0 

GO  TO  (100,200,300,400),  K 
100  CONTINUE 

func=ta*ta+tb*tb 

RETURN 

200  CONTINUE 

C(1  )  =  X(  l)  +  X(2)-2.0 
C(2)  =  X{  1  )*X( 1)-X(2) 

RETURN 

300  CONTINUE 

GC(  1,  1)=1.0 
GC(  l,2)=l.O 
GC<  2,  1)  =  2.0*X(  1 ) 

GC( 2,2)=-1.0 
400  CONTINUE 

GF(  1)  =  2.0*(X(  l)-2.0) 

GE(2)  =  2.0*(  X(?)-1.0) 

RETURN 

END 


SUBROUTINE  FESOIR  ( F , XO , X, N,FBOUND ♦GROUND I 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

10 

15 

c 

c 


20 

21 

C 

C 

25 


26 


REFERENCE.  'NONLINEAR  PROGRAMMING  WITH  THE  AID  OF  A  MULTIPLE- 
GPADTENT  SUMMATION  TECHNIQUE*  BY  W. R.  KLINGMAN  AND  D.M,  HIMMELBLAU 
IN  THE  JOURNAL  OF  THE  ASSOCIATION  FOR  COMPUTING  MACHINERY,  VOL. I  I, 
NO.  A  (OCTOBER,  19641,  PP.  400-415. 

DEFINITION  OF  THE  VARIABLES. 

XO  INITIAL  FEASIBLE  POINT 

X  WORKING  POINT  AND  FINAL  MINIMUM 

S  DIRECTION  VECTOR 

GF  GRADIENT  OF  THE  (UNCTION 

C  THF  CONSTRAINTS 

GC  GRADIENTS  OF  THE  CONSTRAINTS 

DELG  NORMALIZED  GRADIENT  OF  THE  VIOLATED  CONSTRAINT 
DELF  NORMALIZED  GRADIENT  OF  THE  FUNCTION 

NUM  NUMBER  OF  ITERATIONS  THROUGH  POWEL 

N  NUMBER  OF  VARIABLES 

IC  number  OF  CONSTRAINTS 

STEP  VALUE  RETURNED  FROM  POWEL  THAT  GIVES  THE  MINIMUM  IN  THE 
S  DIRECTION 

IFLAG  NUMBER  OF  VIOLATED  CONSTRAINTS 

FHOUNO  CRITERION  TO  SEE  IF  FUNCTION  VALUES  ARE  CONVERGING 
GBOUND  CRITERION  TO  SEE  IF  DIRECTION  VECTOR  IS  CONVERGING 
CODE  CODE  THAT  DETERMINES  WHAT  WILL  BE  PRINTED 

DIMENSION  S( 5) ,GF( 5) ,C( 16» ,GC( 16,51 ,X0(5) ,XI5) ,GSUM(5) ,P(5) 
DIMENSION  0ELG(5) ,DELF(5» 

INTEGER  OBJFN,CQOE 
COMMON/ZANG/  F , I POWEL , OB JEN 
COMMON/CON/IC,KQ,COOE,ICON( 16) 

NUM=0 

START  WITH  FEASIBLE  POINT 
DO  1C  J=1,N 
X( J  »  =  X0( J ) 

VAL  =  FUNC  (4,X,C,GF,GC) 

F  =  FUNC  ( 1 ,X,C,GF,GC) 

SET  DIRECTION  TO  THE  NORMALIZED  GRADIENT  OF  F 
SUM=0. 

DO  20  J=1,N 

FSQD=GF( J)*GF( J) 

SUM*SUM*FSQD 
DO  21  J=l,N 
S(  J )  =  -GF( J )/SORT( SUM) 

COMPUTE  MINIMUM  ALONG  THE  S  DIRECTION 
NUM=NUM*1 

CALL  POWEL  (N, X,S, STEP, F, IFLAG) 

DO  26  J=1,N 

X( J )=X( J )+STEP*S( J) 

IF  (CODt.NE.O)  WRITF(KQ,50C)  NUM , F , STE P , ( J , X ( J ) , J= 1 , N) 
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27  IF  (NUM.EQ.U  GO  TO  28 
C 

C  CHECK  FUR  CONVERGENCE 

IF  (  ABS(F-FSAVE) .LT.FBOUNDI  GO  TO  110 

28  FSAVF=F 

IF  (  IFLAG.EO.CI  GO  TO  15 
VAL  =  FUNC  (3|X,C,GF,GC) 

C 

C  COMPUTE  SUM  OF  NORMALIZED  GRADIENTS  OF  VIOLATED  CONSTRAINTS 

DO  30  J=l,N 
30  GSUM(J)=0.C 

DO  60  I=ltIC 

IF  (  ICON! I ) .EO. 1 )  GO  TO  60 

SUM=0.0 

DO  40  J=l,N 

GSOn=GC( I , JI*GC( I, J) 

40  SUM=SUM+GSQO 

DO  50  J=1,N 

DELG(JI=  GC(I,J)/  SQRT(SUM) 

IF  ((GC(I,J).LT.O.).AND.(OELG(J).GT.O.» I  DEL G ( J ) =-DELG ( J ) 

50  GSUM(J)=  GSUM(  J»«-DELG(  J) 

IF  (CODE.LE.H  GO  TO  60 
WRITE  (K0,540)  1 

WRITE  (KQ,  550)  ( GC ( T  , J ) t DE « G( J I , J= 1 1 N) 

60  CONTINUE 

C 

C  COMPUTE  THE  NORMALIZED  GRADIENT  OF  F 

SUM^C.O 
DO  70  J=1,N 

FSQO=GF(  J)’'GF(  J  > 

70  SUM=SUM+FSQO 

A3DELF=  SQRT(SUM) 

DO  80  J=1,N 

OeLF(J|=  GF(J)/A60ELF 
C 

C  S  IS  THE  NEW  DIRECTION 

S( JI=  -DELF( J  )  -GSUM( J) 

80  continue 

IF  (COOE.LF.l)  GO  TO  85 
WRITE  {K0,560» 

write  (KQ,570I  (GF(J),DELF(J),S<J) ,J=l»NJ 
85  SUM=0,0 

DO  RO  J=l,N 
SSQD=S( J )*S( J ) 

90  SUM=SUM+SSQO 

DELS=SQRT( SUM) 

C 

C  CHECK  FOP  CONVERGENCE 

IFIDELS.LT. GROUND  )  GO  TO  IIC 
DO  100  J=1,N 

C 

C  S  NORMALIZED 

100  S( J)=S( J  )/OFLS 
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110 

120 


500 

510 

520 

530 


540 

550 

560 

570 


4 


8 

5 


20 


25 


30 


32 


GO  TO  25 
WRITE  (Kq,510» 


DO  120  J=1»N 

WRITE  (K0,520) 
WRITE  (KQ,530) 
RETURN 

FORMAT  (1H0,I6. 
1E15.8,/.8H  STEP 
FORMAT  (IHO, 
FORMAT  (4H  XO(i 
FORMAT  (54H  THE 
lE15.B,/f40H  THE 
139H  the  number 


J,XO( Jl , J,X( J) 

F ,NUM,OBJFN 

21H  RETURNED  FROM  P0WFLLf/»17H  FUNCTION  VALUE  = 
=  ,F15.8,/f(3H  X(,I2,3H)  =fE15.8)) 

13HINITIAL  POI NT f 1 5X » 1 IHFI NAL  POINTt//) 

12, 3H)  =,E15.8,5X,2HX( ,I2,4H)  =  ,E15.8,/) 
MINIMUM  FUNCTION  VALUE  WITHIN  THE  CONSTRAINTS 
number  of  ITERATIONS  OF  POWELL  IS  =16,/, 

OF  FUNCTION  EVALUATIONS  IS  ,16) 


FORMAT  (IIH  CONSTRAINT, 12, /,9H  GRAD lENT , 17X , lOHNORMALl ZED) 
FORMAT  ( 1X,E15.8,10X,E15.8) 

FORMAT  (9H  FUNC T  ION  , / ,9H  GRAOI ENT, 1 7X, lOHNORMALl ZED, 16X,3HSUM) 
FORMAT  (  1X,E15.B,10X,E15.B,10X,E15.R) 

END 


, 


IS 


, 


SUBROUTINE  POWFL 
DIMENSION  S( 20) , 
INTEGER  OBJFN 
COMMON/ZANG/  F, 
STEP=0.0 
NUM  =  0 


(N,P0,S,STEP,FF,IFLAG) 

P0(20) ,  Z{20) ,  F(4) ,  V(4) , 
IPOWEL, OBJFN 


R=1.5 

A=0.0 

Q=.l 

GO  TO  5 

Q=Q*.l 

R=l. 

DO  8  I  =  1 ,  N 

POI  I  )=P0( I )+STFP*S( I ) 

CONTINUE 

FA=FF 

DO  20  I  =  1,N 

Z( I  )=P0(  I ) ^0*S(  I ) 

CALL  BNDY  (CO, IFLAG,Z,N) 

IF  (IFLAG.GE.l)  GO  TO  25 
FQ  =  FUNC  I  1 , Z , 0 , GF , GC ) 
0BJFN=08JFN^1 
IF  (FO.LT.FA)  GO  TO  50 
CONTINUE 

iflagi^iflag 

DO  3^  I  =  1,N 

Z( 1 )  =  P0( I )-0*S( I  ) 

CALL  BNDY  ( CO t I  FLAG  ,  Z  ,N ) 
IF  (  IFLAG.EO.O)  go  TO  32 
IF  I STEP.EO.O, )  GO  TO  4 
RETURN 

FNC  =  FUNC  I  1 ,Z ,0,GF  ,GC I 
UHJFN=0BJFN+1 


C0(30) 
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IF ( FMQ.GE . FA.ANO . I  FLAG! .EO.OI  GO  TO  40 
IF  (FNQ.LT.FA)  GO  TO  35 
IF  (STEP.EQ.O.)  GO  TO  4 
DO  33  I=l»N 

33  Z( ! )  =  P0( T  >  +  0*5(1) 

CALL  BNDY  ( CO, I FLAGtZ tN) 

RETURN 
35  B=-Q 

FR=FNQ 
ISW=0 
GO  TO  6C 
40  STFP=-Q 

FF=FNQ 
RETURN 
50  B=Q 

FB=FO 
ISW  =  1 

60  CF=  1 .0 

SUM=1.0 
70  CF=CF*R 

SUM-SUM+CF 

IF  (  I5W.EQ. 1 )  GO  TO  80 
C=-Q*SUM 
GO  TO  90 
80  C=Q*SUM 

qO  DO  100  1=1, N 

100  Z(  I  )  =  P0( 1 )+C*S( I ) 

CALL  BNDY  ( CO , IFLAG,Z ,N) 

IF  (  IFLAG.EQ.O)  GO  TO  102 

101  STEP=B 
FF=FB 
RETURN 

10?  CONTINUE 

FC  =  FUNC  ( 1 ,Z ,n,GF ,GC ) 

0PJFN=0BJFN+1 
IF  ( FC.GT.F3)  GO  TO  IIC 
A=B 
FA  =  FR 
D  =  C 
FR=FC 
GO  TO  70 
no  STFP  =  A 

FF=FA 

IF  ( STEP.EQ.G.C)  GO  TO  4 

RETURN 

END 

SUBROUTINE  BNDY  ( C , I F L AG , X , N ) 

DIMENSION  C(  16)  ,  X(  5) 

INTEGER  CODE 

COMMON /CUN/ I C,KQ, CODE , ICON! 16) 

IFL An=0 
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VAL  =  rUNC  (  ?f  X,C,GF  ,r.c ) 
h  =  riJNC  (  1  ,  X,C  ,GF  ,''-C  ) 

IF  (Cnf)t:.Nr,3)  GO  TO  iC 

WPITF  (K!J,'50)  f 

WFITf  (K(D,^0)  (K,X(K)  ♦K=ltN) 

WPirr-  (KQfVCI  (  1.  ,CIL  I  tL  =  l  f  IC) 

10  no  40  j=i ,  ic 

IF  (  ABS(C(  JII.LT.  l.f'E-l'i)  C(J)=C, 

IF  ((■  (  J  I  .LF  .0.0  )  0,0  TO  30 

ICON! J)=-l 

IFLAO.=  IFLAG+1 

IF  ICnDL.NF.  3)  0,0  to  40 
V-Pior  (KOf«C)  J 
GO  TO  40 
30  IC0N(J)=1 

40  CONTINUE 

SO  FOK^^AT  ;4H  F  =tF15.P) 

60  FflFMAT  (  3H  X(,T7,3h)  =,FIS.8» 

70  FORMAT  (20  C, I?, CIS. 8) 

RO  FORMAT  (IIH  constraint, 12, 12H  IS  VIOLATEO) 

RETURN 
ENO 


no 


RICOCHET  GRADIENT  (28  SUBROUTINES  AND  DRIVER) 
PROGRAMMER; 

J.  Greenstadt  and  R.  T.  Mertz,  IBM/Adapted  for  Use  at  NELC  by 
D.  C.  McCall,  Decision  and  Control  Technology  Division 

PURPOSE: 

To  find  the  point  (x  j ,  X2,  . . x^j)  at  which  the  objective  function 
f(xj,  X2,  • . Xj^)  takes  on  its  maximum  value,  subject  to  the  constraints 
gk(xi,X2. 

RESTRICTIONS  AND  LIMITATIONS; 

The  program  handles  up  to  50  constraints  and  50  varisbles.  The  first 
partial  derivatives  of  the  function  and  all  constraints  must  be  obtainable. 

LANGUAGE; 

FORTRAN  IV 

COMPUTER  CONFIGURATIONS: 

IBM  360/65 

ENTRY  POINTS: 

Main 

SUBPROGRAMS  AND  WHERE  REFClvENCED: 

Programmer-supplied  Programs 

PROB  called  by  OBFUNC,  OBGRAD,  CONSTR,  CSTRNM 

Library  Subprograms 
ABS 
SQRT 

DEFINITION  OF  VARIABLES: 

The  manual  gives  a  complete  definition  of  the  variables  and  gives  a 
summary  for  each  subroutine. 

INPUT  FORMAT: 

A  subroutine  PROB  and  a  set  of  data  cards  are  needed. 

The  subroutirie  defines  the  objective  function,  constraints,  and 
gradients.  The  form  used  is 

SUBROUTINE  PROB  (NUMX,  X.  KK,  INDX,  VAR.  GRAD,  N,  NC, 
C) 


III 


The  necessary  dirneiuioning  is 

DIMENSION  X(50),  GRAD(50),  NC(200),  C(200) 


DEFINITION  OF  THE  VARIABLES: 

NUMX  Seiial  number  of  X  (unused  by  PROB) 

X  The  coordinates  of  the  point 

KK  =  0  Gives  a  value  of  the  objective  function 

1  Gives  a  value  of  the  KKth  constraint 

INDX  =  0  Gives  a  function  value 

1  Gi/es  a  gradient  value 

VAR  Returns  the  objective  or  constraint  function 

vaiue 


GRAD  Returns  the  N-gradients 

N  Number  of  variables 

NC  Can  be  used  to  read  in  integer  data  (up  to  200) 

C  Can  be  used  to  read  in  real  data  (up  to  200) 

See  sample  subroutine. 


SAMPLE  SUBROUTINE 


P'U'.M  (UiiMX»X/KK»  irnX»Vflf{>OP. 

NLO  iHxbtH.liii  i'  j  i 

.•■Atj  <k:id  r)T-'  '‘f.i.mAij  iTE'^t  >-'y/'i’!]LC 


r(5i0!-i  x(5(l)  »o:^''i)(b()i  fNC 


c  •  •  < 

c 

R 

I  JDX)  210»  liOf  120 


C 
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L  VAL 
VAL: 

[ul  1  ! 

k- _ L 

i 

1  I 

i 

._J 

t _ 

1 _ 

* 

9 

IF(KK)  2(10.1000*  112 

IFcKk-Z)  113.i;f)ii2»  200  '  ! 

j 

These  *.i.J6t  a 

igree  with  th 

e  no.  of  con 

Ltralnts 

■ 

nvmm 

om 

JGOl 
GO  1 

=  KK  1 

lo  (inoi .  I(j02)  r  JGOl 

1 

_  i 

one  label  fdr  each  const 

raint. 

liQnvBn 


CORAniFNT  ] 


DO  121  0-i..^  : 


|QrQBn9t|n|iai 


23U»2nOO»  124 
IF(kK-2)  12rjf  c‘n2'..»  210 


0  (::pii)»'’ 020)  >  JG02 


I  i 


■  This  must  a^ree  with  thel  no*  vari 


i  i  I 

;These  must  deree  with  the  no.  of  conktralnCB 


one  label  for  each  set  of  constraint  sradlants 


II 


CALcl'JUAnori  - 

VALJ2S 

— 

/  <  (  X  ( 1 )  + 1 .  )  *  ( X  ( 1 )  + 1 « )  •<•  V  ( 2  1  <  X  ( 2  )  )  ,  Obicctive  function 


Confstraints 


GO  TO  200 


CALCUl.ATIO.J  -  'GRAUIfTl'T 


*12000  lUF.IiO  /.=  (  X  t  1  )  +  l  .  )  ♦  (  X  (  U  +  1  .  )  tvi.?)  .x  (2) 


n>,i 


t  of  o'olective  function 


0  21  n 


Orndients  of  constraints 


GU  To 


*  indicates  cards  chat  are  clianged  from  nrograr.  co  ncogra 
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The  user-supplied  data  follow  the  required  data.  The  user’s  data  are 
divided  into  seven  classes;  preceding  each  clas  ;  is  a  header  card  with  the 
class  number  in  the  first  column.  All  other  data  cards  must  be  blank  in  the 
first  column  and  are  read  only  through  column  72. 


The  seven  classes  arc: 

0  Descnption  of  th;,  prcblern 

1  Setting  for  the  output  conttol  switches  (See  OUTPUT  FORMAT) 

2  Integer  parameters  used  in  the  algorithm  subroutines 

3  Real  parameters  used  in  the  algorithm  subroutines,  including  the  feasible 
starting  point 

Classes  2  and  3  usually  use  the  same  parameters  from  problem  to  problem 
except  for  the  starting  point. 

The  starting  point  begins  in  the  47th  entry  under  Class  3,  and  is  read  in 
with  an  EI4.8,  3E15.8  format. 

4  Size  of  the  problem.  The  numbe**  of  variables  ki  columns  2-5;  the  number 
of  constraints  in  tlumns  6- 1C. 

5  Integer  Parameters  for  the  PROB  subroutine  -  14,  1 515  format 

u  Real  Parameters  for  the  PROB  subroutine  -  E14.8,  3E15.8  format 

The  read!  -6  in  of  data  is  terminated  by  two  cards  with  a  9  in  column  1 .  IT 
tiij  program  is  to  be  run  more  than  once  with  diffeient  data,  one  card  with 
a  9  in  column  1  is  placed  between  the  data  for  the  different  problems  and 
iwo  9-cards  are  u>  cd  to  terminate  the  program. 


SAMPLE  DATA  DECK 
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OUTPUT  FORMAT: 


The  output  for  the  program  is  controlled  by  control  switches  set 
under  the  first  class  of  data.  Following  the  header  card  are  10  cards  numbered 
consecutively  from  1-10  in  columns  4  and  5.  The  switch  controls  are  in 
columns  6  through  72.  A  “1”  punch  indicates  the  value  will  be  printed;  a 
blank  indicates  the  value  will  be  omitted.  In  the  description  of  each  subroutine 
in  the  manual  the  code  for  printing  the  values  is  given.  For  example,  to  print 
the  final  function  value  from  the  subroutine  MONITR,  we  use  the  code  1 , 44, 

F  (on  page  38).  A  “T  in  column  44  of  the  first  card  will  cause  F  to  be 
printed.  This  code  appears  with  the  printout  to  help  with  the  identification 
of  the  values. 

A  good  choice  for  the  minimal  amount  of  printout  Gust  the  final 
result)  is  I’s  in  columns  11,  16,  17,  35,  38,  44,  45,  48,  51,  54  on  the  first 
card;  the  other  cards  are  blank  in  columns  6-72. 

The  maximum  printout  occurs  when  all  10  cards  are  blank  in  columns 
6  through  72.  For  debugging  purposes  the  values  desired  may  be  chosen  from 
the  manual. 

MATHEMATICAL  METHOD  AND  REFERENCE: 

See  the  manual  Contributed  Program  Library  #360D-1 5.3001  ‘Non- 
Linear  Optimization-Ricochet  Gradient  Method. 

TEST  EXAMPLES: 

The  following  problems  were  tested  with  this  program: 

Problem  1. 

Minimize  f  =  (Xj  -  2)^  +  (X2  -  1)^ 

gi  -  Xi  +  X2  -  2  <  0 

g2  =  x“-X2<0 


Problem  2.  (Circle) 

1 

Maximize  f=  ~  ~ 

(x|  +  D"  x:; 

subject  to  g[  =  -  4  >  0 

g  >  =  16  -  x“  -  x'^  >  0 


Ret.  Klinginan  and  Himrnclblau 


Problem  3. 
Maximize 


f  =  y  +  sin  X 


subject  to  0  <  X  <  1 
0<y  <  1 

x^  +  y^  <  1 


W.  -  1-366 


Ref.  problems  3  through  6  were  from  Krolak  and  Cooper  as  found  in  the 
Klingman  and  Himmelblau  paper. 6^ 


Problem  4. 
Maximize 


f=-(y  -x)4  +  (l  -X) 


subject  to  0.2  <  X  <  2.0 
0.2  <y  <2.0 

X^  +  y2  <  I 


fmax.  =  0  8 


Problem  5. 
Maximize 
subject  to 


f  =  -x^  +  x-y3  +  y  +  4 

0.2<x<2.0 

0.2  <y  <2.0 

X^  +  y2  <4 

fmax.  =  4.5 


Problem  6. 
NK'';;:nize 


f  =  exp\-(x  -  1) 


(y3- 0.5)3 

0.132  , 


subject  to  0.2  <  X  <  2.0 
0.2  <y  <2.0 

*1  T 

X-  +  y-  <4 
Imax.  =  1 .0 


Problem  7.  (Lootsina) 


Minimize 


f  =  x'l  -  bx~  +  1  lX[  +  x^ 


RESULTS  OF  TESTED  EXAMPLES; 


1.  Problem  1 

Initial  Point 

Computed 

Optimum 

Objective 

Function 

Number  of 

Iterations 

’‘I 

-1.0 

0.99999994 

-1.000 

47 

^2 

2.0 

0.99999988 

-1.0 

0.99999994 

-1.00 

36 

^2 

0.0 

0.99999988 

^1 

2.0 

0.99999994 

-1.000 

7 

^^2 

2.0 

0.99999988 

2.  Circle 

^1 

2.00 

-2.0 

1.000 

265 

Xt 

im 

-2.8 

0.00076647 

3.0 

2.0 

0.111111 

12 

0.0 

0.0 

5.0 

-1.9999857 

0.999971 

304 

^2 

5.0 

-0.00756613 

’‘I 

0.0 

-1.99999914 

0.999982 

196 

^2 

3.0 

0.0060059 

’‘I 

1.0 

-2.0 

1.000 

316 

^2 

1.0 

-0.000674 
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3.  Krolak  and  Cooper  #1 


Initial  Point 

Computed 

Optimum 

Object've 

Function 

Number  of 
Iterations 

^-1 

0.0 

0.629556 

1.3657408 

33 

Xt 

0.9 

0.776955 

-0,3 

0.63038069 

1.3657379 

29 

X2 

-0.3 

0.77628624 

’‘I 

0.0 

0.63037407 

1.3657379 

30 

X2 

-2.0 

0.77629149 

1.0 

0.62933254 

1.3657408 

53 

X2 

-2.0 

0.77713621 

1.0 

0.63037902 

1.3657379 

30 

X2 

1.0 

0.77628750 

4. 

Krolak  and  Cooper  #2 

’"1, 

-.5 

0.19999999 

0.80000001 

7 

X2 

-.5 

0.19999999 

X 

2.0 

0.19999999 

0.80000001 

25 

X2 

0.0 

0.19999999 

1.0 

0.19999999 

0.80000001 

29 

X-) 

1.0 

0.19999999 

5.  Krolak  and  Cooper  #3 


Initial  Point 

Computed 

Optimum 

Objective 

Function 

Number  of 
Iterations 

0.0 

0.50000012 

4.5000 

36 

^2 

0.9 

0.50000018 

’‘I 

-0.3 

0.50000030 

4.5000 

36 

^2 

-0.3 

0.50000030 

0.0 

0.50000006 

4.5000 

38 

^^2 

-2.0 

0.50000006 

1.0 

0.50000000 

4.5000 

6 

^2 

0.0 

0.49999994 

1.0 

0.50000006 

4.5000 

37 

^2 

1.0 

0.50000006 

6. 

Krolak  and  Cooper  #4 

-0.5 

0.99974561 

0.99999994 

106 

^‘2 

-0.5 

0.70710695 

’‘I 

2.0 

1. 0002337 

0.99999994 

82 

^2 

0.0 

0.70710659 

’‘I 

1.0 

1.0000000 

1.000000 

13 

^^2 

1.0 

0.7071 1035 

7. 

Lootsma 

’‘I 

0.37896395 

0.0 

Xt 

i. 

1.6807594 

1.4142141 

1.4142141 

52 

^3 

2.3471994 

1.4142141 

8.  Colville  1 


Initial  Point 

Computed 

Optimum 

Objective 

Function 

Number  of 
Iterations 

o 

d 

X 

0.30000037 

-32.343491 

337 

X2  0.0 

0.33695012 

X3  0.0 

0.40000027 

X4  0.0 

0.43677974 

X 

b 

0.21579641 
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